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Abstract 

We derive solutions of the Kirchhoff equations for a knot tied on an infinitely long 
elastic rod subjected to combined tension and twist. We consider the case of simple (trefoil) 
and double (cinquefoil) knots; other knot topologies can be investigated similarly. The 
rod model is based on Hookean elasticity but is geometrically non-linear. The problem is 
formulated as a non-linear self-contact problem with unknown contact regions. It is solved 
by means of matched asymptotic expansions in the limit of a loose knot. Without any a 
priori eissumption, we derive the topology of the contact set, which consists of an interval 
of contact flanked by two isolated points of contacts. We study the influence of the applied 
twist on the equilibrium. 



1 Introduction 



Knots are found in everyday life, shoe lacing being probably the most common example. They 
are also essential in a number of activities such as climbing and sailing. In science, knots have 
long been studied in the field of mathematics, the main motivation being to propose a topological 
classification of the various knot types, see the review by Tabor and Klapper (19941. Recently, 



there has been an upsurge of interest in knots in the biological context: knots form spontaneously 
in many long polymers chains such as DNA (Katritch et al. 19961 or proteins, and have been 



tied on biological filaments (Arai et al. 19991. Knotted filaments have a lower resistance to 



tension than unknotted ones and break preferably at the knot (Saitta et al. 1999 Pieranski 



et al. 2001a I. Despite a wide range of potential applications, the mechanics of knots is little 



advanced. The present paper is an attempt to approach knots from a mechanical perspective by 
using a well-established model of thin elastic rods. 

The problem of finding so-called ideal knot shapes has received much attention in the past 
decade (Katritch et al. 1996 Stasiak et al. 19981. In this geometrical description of tight knots. 



a impenetrable tube with constant radius is drawn around an inextensible curve in Euclidean 
space and one seeks, for each knot type, the configurations of the curve such that the radius of 
the tube is maximum. The case of open knots, where the curve does not close upon itself, has 



under tension (Pieranski et al. 2001a I. 



been studied by Pieranski et al. (2001b I in connection with the breakage of knotted filaments 



To go beyond a purely geometrical description of knots, it is natural to formulate the problem 
in the framework of the theory of elasticity. The case of tight knots, or even of moderately tight 
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knots, leads to a problem of 3D elasticity with geometrical nonlinearities (finite rotations), finite 
strains, and self-contact along an unknown surface: there is no hope to derive analytical solutions. 
Numerical solution of this problem raises considerable difficulties too, which have not yet been 
tackled to the best of our knowledge. In the present paper, we study the limit of loose knots, 
when the total contour length captured in the knot is much larger than the radius of the filament. 
In this limit, it is possible to use a Cosserat type model and describe the rod as an inextensible 
curve embedded with a material frame, obeying Kirchhoff equations; as we show, the equilibria 
of open knots can be solved analytically in this limit. 

Self-contact in continuum mechanics, and in the theory of elastic rods in particular, leads 
to problems that are both interesting and difficult. This comes from the fact that the set of 
points in contact is not known in advance — in fact, not even the topology of this set is known. 



This paper builds up on prior work by von der Mosel (19991; Schuricht and von der Mosel 



(2003), who characterizes the smoothness of the contact force in equilibria of elastic rods, and 
by Coleman and Swigon ( 2000 1 , who write down the Kirchhoff equations for rods in self-contact 
explicitly, including the unknown contact force. These equations have been solved by numerical 
continuation in specific geometries by Coleman and Swigon (2000); van der Heijden et al. (2003); 



Neukirch (2004). In these papers, the authors simultaneously solve for the non-linear Kirchhoff 
equations and for the unknown contact forces. In the present paper, we shown that, under the 
same set of assumptions that warrant applicability of the Kirchhoff equations, one can in fact 
neglect the geometrical nonlinearities in the region of self-contact. As a result, nonlinearities 
and contact can be addressed in well separated spatial domains. This brings in an important 
simplification and, as the result, we are able for the first time to derive analytical solutions of a 
self-contact problem for rods undergoing finite displacement, exhibiting a non-trivial contact set 
topology. 

Our solution is constructed with matched asymptotic expansions with respect to a small 
parameter e which is zero for a perfectly thin rod. As is done routinely in boundary layer analysis, 
we use qualitative reasonings (dimensional analysis) to justify how the various quantities scale 
with the small parameter e. We emphasize that our final solution is exact and does not involve 
any other assumption than the smallness of the parameter e : it is asymptotically exact. Our 
presentation is based on formal expansions; proofs of convergence are beyond the scope of the 
present paper and can hopefully be established in the future. 

In the present paper, we consider a knot loaded under mixed tension T and twist U . In a 



previous short paper (Audoly et al. 2007), we have announced some of the results reported here, 
for the case of a purely tensile loading, [7 = 0. In addition to presenting a justification of these 
results, we address here the influence of twist on the knot shape. 

The outline of the present paper is as follows. In Section |2] we introduce the Kirchhoff 
equations for rods in equilibrium, including the contact forces relevant for the knotted geometry; 
we discuss the equivalent formulation as a minimization problem with topological constraints. 
In Section |3] we consider the case of an elastic curve with vanishing thickness and show that 
the region of contact collapses to a point connecting a circular loop and two straight tails. In 
Section|4] we carry out the dimensional analysis of the solution with small but nonzero thickness. 
We show that the equilibrium solution is composed of three types of regions, namely a loop 
and two tails connected by a braid. The scaling of the unknowns with the small dimensionless 
parameter e are identified. Following the general methodology of matched asymptotic expansion, 
we use these scalings to devise a perturbation scheme of the original equations in (non-integer) 
powers of e. The resulting equations are written down and solved in the various regions: the 
tails are solved in Section [5] the loop in Section [6] The solution in the braid region is the 
most challenging as this is where contact occurs, and in Section [7] we obtain a universal solution 
describing the shape of the rod in this region. In Section [8] we build a global solution by matching 
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Figure 1: Two knot types are considered here: (a) simple open knot, also known as trefoil knot, 
noted 3i, and (b) double open knot, also known as cinquefoil knot, noted 5i. The theory can be 
extended to other knot types. 

the solutions derived previously in each region. We obtain a unique equilibrium solution for any 
given value of the loading parameters (force and twist). In Section |9] this theory is validated by 
experiments. Appendix [A| discusses the topology of the contact set in more details. 

2 Model 

We seek equilibrium solutions of a thin elastic rod bent into an open knot with a prescribed 
type, and subjected to tensile end force and torsional end moment, as shown on Fig. [2] In the 
present paper, we focus on two specific knot types, which are open trefoil knots, also called simple 
knot and noted 3i, and open cinquefoil knots, also called double knot and noted 5i, see Fig. [T] 
Other knot types can be handled similarly. The rod is infinitely long and the loading is applied 
at infinity. 

Our model is based on the Kirchhoff equations for the mechanical equilibrium of elastic rods. 
We consider the case of an unshearable |^ inextensible rod with circular cross-section — this is 
the standard model for elastic rods, which can be derived under fairly general hypotheses from 3D 
elasticity theory]^ Contact of the rod with itself is assumed to be frictionless. The mathematical 
formulation of the problem is based on classical models and is relatively straightforward; the 
challenge of the present analysis is to deal with geometrical nonlinearities and self-contact — 
one of our contributions is to determine the topology of the contact set which is not known in 
advance. 

In the present section, we recall the Kirchhoff equations for rods and show how they can be 
applied to the geometry considered. We emphasize the minimization problem underlying the 
equations of equilibrium, and put the equations in a dimensionless form. 

2.1 Kinematics 

We consider an infinite isotropic elastic rod, bent into an open knot as shown in Fig. |2] with a 
circular cross section of radius h, a bending modulus B and a twisting modulus C. Centerline 
of the rod is parameterized by the arc-length s and is defined by its Cartesian equation, 

r(s) = {x{s),y{s),z{s)) , 

^In topology, a knot is defined as a closed, non self-intersecting curve. Here we consider curves having two 
infinite tails, hence the name 'open knots'. 

^By unshearable, we mean that the rod satisfies the Navier-BernouUi kinematical hypothesis. 
■^Extensions of the present results to different rod models do not raise any fundamental difficulty. 
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where the orientation of the axes is specified below. The tangent to the centerhne is noted 

t(.) ^ % (1) 

Since the rod is assumed inextensible, the tangent is a unit vector, 

|t(s)| = l (2) 

for all s. We note m(s) the internal moment in the rod and n(s) the internal force — these 
variables describe stress distribution in the cross-section in Kirchhoff theory of rods. 



2.2 Constitutive relations 

We assume a linear elastic response (Hookean elasticity) , which is consistent with the small strain 
approximation underlying Kirchhoff theory. The constitutive law for a rod with symmetric (e. g. 
circular) cross-section can be conveniently written in vector form ( jLandau and Lifshitz 1981[ ): 



m(s) = Bt(s) xt'(s)-|-CT(s)t(s), (3) 

where t(s) is the material twist of the rod. The above expression is a condensed form of the 
constitutive relations for a rod that are usually written in coordinates in the material frame. 
The first term in the right-hand side is the bending moment and lies in the cross-section; for 
a symmetric rod, this bending moment is the binormal, t x t', times the bending stiffness B. 
The second term in the right-hand side is the twisting moment and is along the tangent: the 
twisting moment is the material twist, r, times the twist stiffness C. Since the rod is considered 
inextensible and unshearable, the internal force is the Lagrange multiplier associated with these 
kinematical constraints, and it not given by a constitutive law. 



2.3 Loading 

At the end of the rod corresponding to s ^ +oo, a tensile force T and a torsional moment U 
are applied, see Fig. [2] These two vectors are assumed to be coUinear, and are used to define the 
axis z. Global mechanical equilibrium requires that an opposite force — T and moment — U are 
applied at the other end, s — > —oo. At equilibrium, the two long tails of the rod will be aligned 
with the direction z of the force. Owing to our choice of axis, we write 

T = Te2 and \i = U a,^. 



Stability of the long tails require T > but the twist U can be positive or negative. 



2.4 Symmetry 

Given the symmetry of the loading, we focus |^ on equilibrium solutions that are symmetric. 
More accurately, we assume that the knotted rod is invariant by rotation with angle tt about an 
axis perpendicular to the axis z defined by the loading. This is consistent as the endpoints at 

*Since the equations are non-linear, one could argue that some solutions having no symmetry at all could exist, 
as happens in buckling problems. We would miss such solutions since we restrict the analysis to the symmetric 
case from the beginning. This remark applies to the case of a finite thickness h, but to the problem addressed 
here, which concerns the limit of a small thickness h. In this limit, we show that the equations can be linearized; 
using the principle of linear superposition, one can focus on symmetric solutions. 
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Figure 2: An infinitely long rod is bent into a knot with a given type, here a trefoil knot (3i), and 
loaded with combined twist U and axial force T. In this paper, we derive equilibrium solutions 
for this non-linear self-contact problem. 



infinity are swapped by this transformation, and so the loading is globally invariant under this 
transformation 

Let us call y the axis defining this symmetry by rotation with an angle n. The intersection of 
the perpendicular axes z and y defined so far will be the origin O of our Cartesian coordinates. 
The direction perpendicular to y and z defines the third axis x, in such a way that (x,y,z) is 
direct and orthonormal. Intersection of the axis of symmetry y with the centerline defines what 
can be called the midpoint of the rod — intuitively, this is the bottom of the loop in Fig. |2] 
This midpoint is taken as the origin of the arc-length coordinate, s = 0. With this convention, 
the symmetry by rotation about y with angle tt maps a point on the centerline with coordinate 
s onto the point with opposite coordinate (— s). Using this property, it is sufficient to find the 
equilibrium shape of the rod over one half, say the positive half < s < -l-oo: the other half can 
be found by applying the symmetry. 



2.5 Variational formulation, constraints 

The equilibrium shape of the knotted rod can be found be solving a minimization problem: this 
equilibrium shape is a minimizer of the total energy of the rod (potential energy associated with 
loading at endpoints plus elastic energy) under the combined constraints of inextensibility, non- 
penetration and prescribed knot topology. This variational view of the problem will be useful 
later for solving the braid region in Section [7] 
The total energy of the rod is defined as: 

(-K^ + -T^\ds + TD^-URo^, (4) 

where k and r stand for the curvature and the twist of the rod. The integral term in the 
right-hand side is the elastic energy associated with the constitutive law The bending term 
depends on the scalar curvature 

^=\t'{s)\. (5) 

^Note that the knot is not invariant by a reflection with respect to a plane perpendicular to the z axis: this 
reflection leaves the loading globally invariant but changes the knot type, turning a left-handed knot into a 
right-handed one. 
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The last two terms in equation Q represent the work of the apphed tensile force Te^ , related to 
the end-to-end shortening Dao, and of the applied torsional moment Ubz, related to the relative 
rotation i?oo of the ends. 

The minimization of this energy is subjected to a series of constraints. First, the inexten- 
sibility constraint is expressed by equation Q . Second, the topology of the knot is prescribed 
(this topological constraint cannot be written down easily in the general case; it will be shown to 
impose the value of a winding index in the braid when we focus on loose knots later on). Third 
and lastly, one has to consider the non-penetration constraint which can be expressed as: 

|r(si)-r(s2)| >2/i, (6) 

for any si and S2 such that |si — S2I > 4/i. Note that the radius of the rod h enters in the equation 
at this point in the right-hand side of equation The trick of restricting the penetration test 
to couples of points (si,S2) separated by a curvilinear distance greater than 4/i is due to 



von 



der Mosel (19991, and avoids mistaking close neighbors on the centerline for points violating 
the non-penetration condition — it is given for mathematical consistency but is not needed in 
the following: for the problem we consider, we know a priori that the arc-length separation of 
two points in contacts is large, namely of order 27r R where the radius R of the loop is a known 
quantity of order 1. 



2.6 Equilibrium: Kirchhoff equations 



The equilibrium equations for a rod can be derived from the energy ^ by the Euler-Lagrange 
method (Bourgat et al. 1988 Steigmann and Faulkner 19931. This leads to the following 
equations: 



r'(.)=t(.) 

, m(s) , , 
t'is) = X t(,s) 

m'(s) +t(s) X n(s) = 
n'(s) +p(s) = 



(7a) 
(7b) 

(7c) 
(7d) 



where primes denote derivation with respect to arc-length s. The first equation is the defini- 
tion of the tangent, already encountered in equation Q. The second equation combines the 
constitutive relations (|3| with the definition ([5| of curvature. The last two equations express 
the equilibrium of moments and forces on an infinitesimal rod element, and are known as the 
Kirchhoff equations ( [Landau and Lifshitz [l98l] ) . The vector p(s) is the density of distributed 
force per unit length applied on the rod, sometimes referred to as the contact pressure. In the 
present problem, gravity is neglected and the only force p(s) to be considered is the one arising 
from the contact pressure in the regions of contact — if there is no contact, p(s) = 0. 



2.7 Contact set, contact force 

Let us define the contact set as the set of couples of arc-lengths, (si, S2), defining cross-sections 
that are in contact: 

€^{{si,S2) such that |si - S2I > 4/i and |r(si) - r(s2)| = 2 /i}, (8) 

where the first inequality, \si — S2I > 4/i, is to avoid mistaking close neighbors for regions of 
penetrations, as explained earlier, and the second inequality |r(si) — r(s2)| = 2 ft, is the contact 
criterion. 
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We are touching here the main challenge of self-contact problems: the profile of the contact 
pressure p(s) is required to compute the centerline by integration of the KirchhofF equations, but 
it depends itself non-linearly on the geometry of contact, that is on the shape of the centerline. 
In other words, p(s) and the contact set £ must be determined in a self-consistent way but none 
is known a priori. In particular, the topology of the contact set is not known in advance. It will 
be obtained later as an outcome of our calculations. 

For any couple (si,S2) in the set £, the corresponding cross-section are in contact. By the 
action-reaction principle, we have p(si) — — p(s2)- In addition we assume frictionless contact: 
the force p has to be normal to the rod surface, and so is aligned with (r(s2) — r(si))- This 
implies that the contact force is aligned with the vector joining the points r(si) and v{s2)'- 

Ms,) ~ vis,)) f{x{sr)-x{s2))/{2h)\ 
p(.i)^p(.i) ^ ^ " ^''» ^p(.i) {y{s^)-y{s,))/{2h) , (9) 

\{z{s^)-z{s,))/{2h)J 

for (si,S2) G C In this equation we have introduced the scalar contact pressure p{s); since the 
rod is a ID object, the contact force has the dimension of a force per unit length but we shall 
nevertheless call it a contact pressure. For the solution to be physical, the pressure must be 
positive: 

p{s) > 0. (10) 
In terms of the scalar contact pressure, the action-reaction principle can be rewritten as 

p{si)^pis2). (11) 

2.8 Boundary conditions 



Thanks to the symmetry introduced in Section 2.4 the equations of equilibrium (|7| need be 



solved over half the rod only, that is for < s < +oo. These equations form a boundary 
value problem as there are conditions to be satisfied at both ends of the interval. The following 
conditions must be satisfied at the endpoint s = +oo: 

m{+oo) = Ue^, (12a) 

n(-l-oo) =Te^, (12b) 

t(+(^)=e„ (12c) 

r{+oo) X t(+oo) = 0. (12d) 

The asymptotic conditions holding at the opposite end of the rod, s — > —oo, can be found by 
symmetry. The first two asymptotic conditions above enforce the loading applied at infinity. 



Condition ( 12c I imposes that the tangent is asymptotically aligned with the applied force, which 



is an obvious necessary condition for minimizing the energy. The last condition ( 12d I defines 
the z axis as the asymptote of the centerline far away from the knot — without this convention 
there would be infinitely many solutions, corresponding to the invariance of the system under 
rigid-body translations perpendicular to the z axis. 

The boundary conditions at the midpoint s = of the rod derive from the invariance of the 
solution by a rotation of angle tt about the y axis: 

(13a) 
(13b) 
(13c) 
(13d) 



t(0). 


By 


= 


m(0) • 


By 


= 


n(0). 


Gy 


= 


r(0) X 


By 


= 
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The justification for each of these three equations is similar, and will be given here for the first one 
only. Let us write the Cartesian coordinates of the tangent t(0) at midpoint as t(0) = (^§,^0,^0)- 
According to our symmetry assumption, rotating the system by an angle tt about the y axis is 
equivalent to reversing the orientation of the centerline. The first operation changes the tangent 
to (— tg, +iQ, — tp), while the second one changes it to its opposite, (— tg, — ig, — tg). Equality of 



these two vectors imposes = 0, which yields equation (13a I 



2.9 Invariants 

Due to their variational nature, the equilibrium equations ([t]) are associated with several in- 
variants as discussed by Maddocks and Dichmann ( 1994[ ). These invariants are known to be 



conserved in the absence of distributed force, p = 0. In the present case, the distributed force 
can be nonzero but remains everywhere perpendicular to the tangent t{s). Under this assump- 
tion, it is straightforward to check that the following expressions are still invariants: 

/i = m(s).t(s) and = + n(s) • t(s). (14) 

The first invariant is directly proportional to the material twist r(s) — Ii/C, and is known to 
be uniform for a rod with symmetric cross-section in equilibrium in the absence of distributed 
torques. The famous KirchhofF analogy identifies the equations of equilibrium of a symmetric 
rods to the equations of motion of a table top subjected to gravity. According to this analogy, 
the first invariant expresses conservation of the angular moment about the axis of the top; the 
second invariant expresses conservation of the energy of the top. 



The constant values of these invariants are imposed by the boundary conditions (12 1 



Ii = U = Ct{s) and therefore t(s) — ^, (15a) 

and 

^2=^ + T. (15b) 
Note that conservation of these invariants requires frictionless contact. 

2.10 Dimensionless form 

The problem has been formulated so far in terms of the loading parameters T and U, of the 
thickness h, and of the elastic stiffnesses B and C. In this section, we use dimensional analysis 
and rewrite the equations in a form that depends on two dimensionless parameters, U and e, 
only. 

To this end, we introduce the characteristic length L* , force F* and moment M* as follows: 

L* = F* = T, M* = VbT. (16) 

These quantities are used to define dimensionless variables, noted with a bar. For instance, we 
define the dimensionless arc-length s and position r as 

f(s)= ^^=r(s)-^-. (17) 



8 



Note that rescaled functions, such as r, are always considered to be a function of a rescaled 
argument, here s and not s: a prime on a barred function implies that derivation is with respect 
to the rescaled arc-length. For instance, the rescaled tangent is defined by 

and happens to be the same unit vector as the physical tangent t, evaluated at the corresponding 
point s = s L*. Similarly, the rescaled curvature is defined as k = |t |: 

k{sL*) 



k{s) = 



1/L* 



The rescaled internal moment and torsional couple are: 

_ mjsL*) mjsL*) U 

mis) = = — , and U = 

^ ^ M* JWT M* 



U 



'BT 



(18) 



(19) 



The internal force n is naturally rescaled using the typical force F* ~ T, while the contact 
force per unit length, p, is rescaled using the dimension F*/L*: 



njsL*) pjsL*) B 

n(s) = — — and p(s) = — - — -. 



(20) 



Having defined the rescaled form of the various quantities, we proceed to rewrite the equations 
of the problem in dimcnsionlcss form. We start with the constitutive relation 



m(s) = t(s) X t (s) + Ut{s). 

The kinematical relations and equilibrium equations Q write; 

r'(s)=t(s), 
t'(s) — m(s) X t(s), 
m'(s) + t(s) X n(s) = 0, 
n'(s) +p(s) = 0. 



For the asymptotic conditions ( 12 1 we obtain: 



m(-l-oo) = U Gz, 
n(+oo) = e^, 

t(+CX)) = Gz, 

r(+cx)) X t(+oo) = 0, 



(21) 



(22a) 

(22b) 
(22c) 
(22d) 



(23a) 
(23b) 
(23c) 
(23d) 



while the midpoint conditions take the same form as the original expressions (13 1, with original 
variables replaced by barred ones: 



t(0) • By = 0, 

m(0) ■ By = Q, 
n(0) -6^ = 0. 
r(0) X = 



(24a) 
(24b) 
(24c) 
(24d) 
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In rescaled form, the invariants (15 1 read: 



Ii = m(s) • t(s) 
7 m2(s) 

^9. — Z 



+ n(s) • t(s) 



U 
~2 



(25a) 
(25b) 



as 



We now turn to the non-penetration constraint. In terms of r = v / L* , Eq. (|6| can be rewritten 

h 



|r(Si)-r(s2)|>2 



L* 



In the right-hand side a fundamental dimensionless number of the problem has appeared, namely 
the ratio of the rod thickness to the characteristic length built from the traction force T and 
the bending stiffness B of the rod. It will be convenient to deal with this dimensionless number 
using an auxiliary number e defined as 



= 21/4 



2/i2r 

B 



1/4 



(26) 



This details of the present definition of e will be motivated later in Eq. (31 1. In terms of e, the 
non-penetration condition can be written as 



(r(si)-r(s2))'>2. 



(27) 



In equations (21 1 to (27 1, we have rewritten all the equations of the problem in terms of two 
dimensionless parameters only, U and e, defined in Eqs. (19) and (26 1 respectively. The first 
parameter U is the rescaled torsional moment; the second parameter e is the aspect-ratio of the 
rod. In this paper, we focus on the limit of thin rod, or a loose knot, e — > 0, and build an 
asymptotic solution of the set of equations above, for axbitraxy values of U . 



3 Zero radius solution 

Before we build a solution of the equations for small e, it is useful to consider the limit of an 
infinitely thin rod, /i = 0, that is the problem of tying a knot on an elastic curve. This case 
corresponds to e = 0, and is the subject of the present section. 

As we shall show later, the limit e — > is singular: the contact region has nonzero length 
for e > but shrinks to a point in the limit e = 0. It follows that the solution of the limit 
problem (e = 0) is less regular than the solutions for e > and must be defined in a weak sense. 
In addition, the topological constraint on the knot type is not easy to write down for the limit 
problem. We shall not study convergence for e — > in a mathematically rigorous way. Following 
a pragmatic approach, a solution for e = is constructed in the present Section based on some 
assumptions (planarity, existence of a point-like contact). These assumptions will be validated 
later when we show that it is possible to extend this e = solution into a family of smooth 
solutions indexed by e > 0. 

3.1 Explicit solution 

In the limit of zero thickness /i = 0, we consider a solution made up of a circular loop connected 
to two straight, semi-infinite tails. Owing to the assumed symmetry of the solution, we focus 
on one half the rod and consider a half-circle starting from the midpoint, connected to a single 
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Figure 3: Case of zero thickness, h — 0. The equihbrium solution is planar and made up of two 
semi-infinite straight tails connected to a perfectly circular loop with radius R. The top of the 
loop is connected to the tails across a singular point O, shown in gray, where both the internal 
force and moment are discontinuous. 



straight, semi-infinite tail. There is no contact, except at the singular point O where the loop 
and tail merge, see Fig. |3] By our previous definition of axes, this point O is the origin of the 
Cartesian frame, the loop is contained in the (?/, z) plane and the tails lie along the z axis. 

Any quantity pertaining to the limit of a zero thickness, h — or e — 0, introduced in the 
present Section will be denoted with a superscript '0'. A subscript 'L' refers to the loop region. 



and T' to the tail region. Let R be the radius of the loop, which will be given in Eq. (30 1, and 
R = R/L* the rescaled radius. The loop region is given by the classical circular solution of the 
Kirchhoff equations. The centerline is a circle given in parametric equation: 

fi(s) = {0, - (R + R cos{s/R)) , ~R sin(s/i?)) , (28a) 
t° (s) = (0, sin(s/i?), - cos(s/i?)) . (28b) 

Note that the constants of integration are such that the top of the loop, s = tt i? is at the origin: 
r5,(7r_R) = 0. The tangent at midpoint is t^(0) = — e^, and that at the top is t"(7ri?) = +6^. 
The internal force and moment in the loop can be found by plugging these expressions into 



equations (21 24 1 with p = 



m5,(s) = (1/i?, U sin(s/i?), cos(s/i?)) , (28c) 
n° (s) = (C7/i?, 0, 0). (28d) 

These solutions describe the loop region, — tt R < s < tt R. 

The solution in the tail is even simpler and describes a straight rod under combined axial 
tension and twist: 

(29a) 
(29b) 
(29c) 
(29d) 

these expressions being applicable for s > tt R. 





= (s - 7TR)e^, 


4(5) 


= ez, 


mi^(s) 


= Uez, 




= ez, 
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3.2 Singular braid point 



At the point connecting the loop and tail regions, s — n R, the solution is discontinuous. Across 
this point, the internal force jumps from n^(7ri?) = U e^/R to ny(7ri?) = e^, while the bending 
moment, defined as the cross-sectional projection of m, drops from e^/R to 0. These discontinu- 
ities point to the presence of contact forces in this region, and will be explained by our analysis 
of the braid region for finite e, see Section [7] 

For this solution to be complete, there remains to compute the radius R of the loop. This 
can be done by writing the conservation of the second invariant given by Eq. ( 25b I across the 
singular point: 



u- = 



u 

~2 



1, 



where the left-hand side comes from the loop region and the right-hand side from the tail. This 
implies 

thatisi?=yj. (30) 



This result was previously obtained by Arai et al. (19991 based on energy minimization of the 
energy (|4| with respect to R. 



4 Perturbation scheme 



In Section |2.10| the equilibrium of a knotted rod has been written as a system of coupled, non- 
linear, ordinary differential equations depending on two dimensionless parameters, \J and e. In 
this paper, we consider the limit of a small e, e <C 1. This limit is in fact the only one consistent 
with Kirchhoff (or Cosserat) description of the rod as a ID elastic object. Indeed, Kirchhoff 
theory comes from a reduction of 3D elasticity, and is justified when the thickness h is much 
smaller than the typical radius of curvature of the centerline. This typical radius of curvature is 
L*, meaning that Kirchhoff approximation makes sense in the limit h <^ L* , that is e <C 1. The 
opposite limit of a perfectly tight knot defines a geometrical problem which has extensively been 
studied, 



Pieranski et al. (2001b I; Katritch et al. (19961; Cantarella et al. (2005) 



The limit e <C 1 under consideration corresponds to a rod whose radius h becomes infinitely 
small while its elastic moduli are kept constant, or equivalently to the case of a fixed radius h 
and elastic moduli when the pulling force becomes very small. We refer to this limit generically 
as the limit of a loose knot. Our somewhat arbitrary definition of the perturbation parameter e 



in Eq. ( 26 1 has in fact been motivated by the simple relation 



(31) 



where h is the rod thickness and R the loop radius defined in Eq. ( 30 1 . The limit of a loose knot 
corresponds to e ^ 0. 

In Section [3) we introduced a solution corresponding to the limit h = 0, that is to e = 0. This 
solution features a singular point where some contact occurs. One of the main contributions of 
the present paper is to come up with a detailed description of this contact region, called the 
braid later on, for small but finite h. A key remark, formulated by [Gallotti and Pierre-Louis] 
(2007), is that the contact region remains very localized for small e. Together with the explicit 
solution for /i = given in Section |3] this suggests the decomposition of the knot solution in three 
domains shown in Fig. |4] a quasi-circular loop, two quasi-rectilinear tails and a braid region in 
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braid 




Figure 4: In the limit of a loose knot considered here, e ^ 1, the equilibrium solution can be 
decomposed into three domains: an almost circular loop, two almost straight tails, and a braid 
region where self-contact takes place. In the vocabulary of asymptotic analysis, the braid region 
is an inner layer, with typical length ^ much smaller than the typical size R of the loop and tail 
regions (outer layers). Note the existence of so-called intermediate region, in darker gray, at the 
overlap between braid and tails, and between loop and braid. 



between. We shall now study the orders of magnitudes of the displacement relevant to these 
different regions. This is an important preliminary step for the quantitative analysis presented 
in the following sections. A simple scaling argument, given in our preliminary paper (Audoly| 



et al. |2007) , shows that the size £ of the contact region is of order of the geometric mean of the 



loop radius R and the rod thickness h, which we write 

VhR. 

This defines an intermediate length scale, between the 'large' length R and the small length h. 
To justify this scaling, we note that the transverse displacement in the braid is fixed by contact 
and is of order h; over a typical length £, this yields a typical curvature h/i'^. At the exit of 
the braid, this curvature has to be matched with that in the loop, of order 1/R. Balancing h/i'^ 
with l/R yields £ ~ VhR as proposed above. 

By this argument the solution features three widely different length scales, R ^ £ ^ h. The 
large scale R is relevant in the loop and tail regions. In the braid region, the relevant length scale 
is £ along the braid axis, and h in the perpendicular direction. Defining the rescaled, typical 
braid length £ and the rescaled radius h by 

we note the orders of magnitude associated with the three fundamental lengths in rescaled form: 

In the vocabulary of inner or boundary layer analysis, the loop and tail regions are both called 
outer regions, while the braid is called the inner region |^ 

^Outer regions are those that are present in the zero thickness solution of Section [s] they undergo a regular 
perturbation for small but nonzero e. In contrast, the inner region (braid) is undefined in the e = solution and 
so has to be built from scratch when e > 0. 
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The above argument clearly shows that the limit e — > is singular and that a uniform 
expansion of the solution with the parameter e is not possible. This is typical of boundary layer 
problem — or inner layer problems in the present case. The classical approach to such problems 
is to use matched asymptotic expansions, that is to build a solution domain by domain using 
different approximations in the outer and inner domains, and to match these solutions in the 
regions of overlap between two adjacent domains. 

As mentioned earlier, the outer regions undergo a regular perturbation. This suggests the 
following, simple expansions in the tail region (subscript T) and in the loop region (subscript L): 




rris) = fUs) + e yris) \ +■■■ and rL(s) = r^(s) + e yi(s) + . . . (32) 




The functions and relevant to the zero thickness case have been given in Eqs. (28a 




and (29a). The six unknown functions xt{s), j/T(s), • • • , zl{s) describe the first-order pertur- 
bation in the tail and loop regions, and will be found later by solving the linearized Kirchhoff 
equations. 

For the inner region (braid, with subscript B), there is no solution available in the limit of 
zero thickness and the above scaling argument suggests an expansion of the form: 

VBia)^ I eTB\ + le^mi^) \ +... (33) 



The first term in the right-hand side represents an infinitesimal rigid-body translation along the 
y axis: the center of symmetry of the braid does not need to remain at the origin when e is 
nonzero, but can only move along the y axis due to the symmetry. The second term is not a 
rigid-body motion. It makes use of the stretched coordinates x — x/e^, y — iij — €.tb) / and 
(7 = z/e suggested by the above scaling analysis: the axial dilation factor 1/e and the transverse 
one comes from the lengths scales I ~ e and and h ^ found earlier. The use of stretched 
variables is classical in problems of elasticity with a small parameter, such as those that arise in 
the analysis of slender elastic bodies. The number tb and the functions xb and ijB are unknowns 
which will be determined later. 

Note that we use the stretched coordinate a as the parameter for the centerline in the braid 
region, and this a = z/e is not the arc-length (in the absence of ambiguity it is common to 
use the same notation for the functions mapping arc- length s to centerline position r, or 
stretched axial variable a to centerline position r). The above scalings imply that the braid is 
almost parallel to the z axis, and so the tangent can nowhere be perpendicular to e^: in the 
braid region, there is a one-to-one mapping between the arc-length and the parameter ct, which 
is proportional to z, and it makes sense to use ct as a parameter along the braid. 



The equations (32 1 and (33 1 provide a starting point for our analysis. These expansions will 
be plugged into the general equations for the knot derived earlier in Section [2] The resulting 
equations for the perturbed tail will be solved in Section [5j those for the perturbed loop in 
Section [6] finally, the leading-order braid solution, which is more difficult to derive, will be 
given Section [7] As implied by the name 'matched asymptotics', the last step is to match the 
solutions obtained in the different regions; this is done in Section [8] by requiring consistency of 
the expansions coming from the two adjacent domains in the regions of overlap. This provides 
a smooth solution of the Kirchhoff equations over the entire domain, which will be shown to be 
unique under some hypotheses. 
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5 Tail solution 



In this section we solve the linearized KirchhofF equations in the tails, which are given by an 
infinitesimal perturbation near the straight solution, see Eq. (32 1. These linearized equations 
are classical and arise in the analysis of linear stability of a straight, twisted rod under helical 
buckling. We characterize the first-order perturbation to the straight configuration due to the 
presence of the knot by computing the functions {xt, Vt, zt). As explained earlier, we focus on 
the tail located on the positive side of the z axis. 



5.1 Linearized Kirchhoff equations near a straight configuration 



To start with, let us plug Eq. (32 1 into the definition ( |7a[ ) of the tangent and compute 

|tT(s)P = |e^|2 + 2ee, •(x^(s),y^(s),z^(s)) + --- = l + 2ez^(s) + ••• 
where the dots stand for higher order terms in e. By the inextensibility condition ([2|, the 
left-hand side has to be equal to 1 for any value of e and so z^(s) = in the right-hand side: 



zt{s) 



(34) 



where is a real constant. This constant can be interpreted as an infinitesimal rigid-body 
translation of the tail along its axis, accommodating the change of curvilinear length captured 
in the loop and braid regions. 

There is no contact in the tail regions and so the contact force p is zero. By Eq. (22d), the 
internal force n is then uniform over the whole tail. Now, this constant value of the internal 



force is set by the asymptotic condition ( 23b I , and so 



(35) 



a quantity that does not depend on e. This makes is possible to integrate the equation for the 
equilibrium of moments: 

mT(s) = + X vt{s), 

where is a constant of integration whose value, ra.K — U ez, is provided by the boundary 
conditions ( 23 1 : 



mr(s) 




(36) 



Equation (22b I for the rate of rotation of the tangent is automatically satisfied at order zero; 
at first order in e, it writes: 

_ ^x'j^isy 

UBzXc \ y'rpis) I +e 

This vector equation is automatically satisfied along the z axis. Projection along x and y axis 
yields a system of two equations: 





x'^{s) - xt{s) + U y'rpis) = 0, 
yris) - i/t(s) - U x'rp{s) = 0. 



(37a) 
(37b) 
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Here we have written the equations for a rod in the small deflection approximation. In the 
problem at hand, it turns out that the tension term dominates over the bending term in the 
balance of transverse forces; this explains why we have a second-order equation rather than the 
classical fourth-order equation of beam problems. 

These equations for a twisted rod linearized near a straight configuration are identical to 



the ones obtained in the linear analysis of helical buckling, see van der Heijden and Thompson 



(2000). Eqs. (37 1 can be put in a compact form when expressed in terms of the complex variable 
wt(s) = xt{s) + i yris): _ 

wt{s)" -WT{s)-iUw'rr{s)^0, (38) 

where i'^ — —1. We seek solutions of this linear differential equation with constant coefficients 
in the form of exponential functions wt(s) = Y e'^'^^^'^ where k and T are complex constants. 
Note that we are free to incorporate the constant term — fc tt i? in the argument of the exponential; 
this amounts to change the definition of the undetermined constant F and will turn out to be 
convenient later on. The possible values of the complex number k are given by the roots of the 



characteristic polynomial of Eq. ( 38 1 : 

k^ -iUk- 1 ^ 0. 
These roots are noted ki = —a + ib and k2 = a + ib where: 




(39) 

The general solution of the equation for wt could be written as a linear combination of the 
functions 

e"°* (cos(6s) + i sin(5s)) and 6+""* (cos(5s) + i sin(6s)) 

but, as said above, it is more convenient to use s — nRas argument. Without no loss of generality, 
we write the general solution as 

wt{s) = r_ e"" t^-'^^) (cos[6 {s-TTR)]+i sm[b {s - ttR)]) 

_^P^g+a(5-^7?) (^cos[b{s-TTR)] + ism[b{s~TTR)]) . (40) 

Here, T_ and r_|_ are two complex constants of integration. The exponentially large solutions 



are incompatible with the boundary conditions (23 1 and so are discarded: = 



Noting A and /i the real and imaginary parts of the unknown complex amplitude r_ = A + i/i, 
we can write the general solution for the displacement as 

xt{s) = n{wT{s)) = (A cos[b{s-TrR)] - n sm[b {s - ttR)]) g-^^^"'^^), 
yrr{s) ^ ^wt{s)) = (m cos[6(s-7ri?)] + Asin[6(s-7ri?)]) e^^^^"'^^). 



By equations ( 29a I , ( 32 1 , and ( |34| the arc- length s is related to the z coordinate by 

z — zt{s) = {s — TT R) + e z!j^ + . . . 

We can use this relation to introduce a change of variable and parameterize the centerline by z 
instead of s. For instance, the above expression for xxis) can be rewritten xt{z) = (A cos(6z) — 

^The value of the tension T 7^ does not appear in the equation as it has been effectively set to 1 by our 
choice of dimensionless variables. 



16 



H sin(6z)) e~°^ + 0{e), where the 0{e) notation means that the equahty is exact up to terms of 
order e. This leads to the following parameterization of the tail, which is valid to first order in e 
included: 

vxiz) = ze^ + e {xt{z) + yriz) ey) -\ 

where 

xt{z) = (A cos(&z) sin(6z)) e^""" H (41a) 

yrp{z) = (/i cos(6z) + A sin(6z)) e""^ H (41b) 

This solution depends on two real parameters, A and fi, which we call the internal parameters of 
the tail. They are referred to collectively as ^t- 

*T = (A,/i), (42) 
and will be determined later by matching with the other regions. 



5.2 Asymptotic expansion near junction with braid 

The matching problem, studied later in Section [8] is based on the expansion of the tail solution 
given above in Eqs. (41 1 near the junction with the braid, that is near the origin z = 0. This 



expansion is computed here. 

The rescaled x and y coordinates of a current point if on the centerline are noted xt and 



yrp. By Eqs. (41 1, their expansion is of the form 

xt{z) = eXT + eX'rrZ + 0{e^,ez'^) 
y^iz) = eYT + eY^z + 0{e^,ez^). 



(43a) 
(43b) 



As implied by the 0(.) notation, the right-hand side is the beginning of an expansion where we 
have neglected terms of order coming from the next order in the global expans ion w ith respect 
to e, and of order ez^ coming from quadratic terms in the expansion of Eqs. (41a) and (41b) 
with respect to z. 



In equations (43) above, the four coefficients {Xt, X^,Yt,Y^) are found by identification 
with the series expansion of xt{z) and yxiz) given in Eqs. (41 1 near z = 0: 





— xt{z 


-0) 


= A 


Yt 


= yriz 


= 0) 


= M 






= 0) 


= —a X — b fi 


Y^ 




= 0) 


= b X — a 11. 



For the matching problem studied later, it is convenient to put these expressions into matrix 
form: 

fXT\ / 1_ 0_\ 

= Mt(U)-^t where Mt(C7) = "^^^^ 



Yt 





V biU) 



1 

-a(U)J 



(44) 



Note that this equation defines the matrix Mt(C/) explicitly as a function of the loading param- 
eter U = U I \/ BT. This matrix 'M.t{U) captures the elastic response of the tail to perturbations 
applied at its end z = 0; it is the only quantity relevant to the tail that will be used in the 
matching problem of Section [8] 
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5.3 Helical instability 



We mentioned that the linearized equations (41 ) arise in the classical analysis of linear stability 
of a straight, twisted rod under helical buckling. For U = ±2, a{U) = in Eq. (39 1 and the two 
complex roots fci and k2 collide: this is the threshold of linear stability for this helical buckling 
mode. This instability will show up in the analysis of the knot later on. 



6 Loop solution 



In this section, we solve the Kirchhoff equation linearized near the planar, circular configuration 



relevant for the loop region. This problem comes from applying perturbation (32 1 to the zero 



thickness solution (28a I for the loop. It is somewhat similar to the classical analysis of stability 



of a circular rod under twist, known as Michell's instability, see Michell (1890|. We focus on 



one half loop, corresponding to the interval < s < tt i?. Deformation of the other half can 
be found by symmetry. The point with arc-length s = is the bottom of the loop and that 
with coordinate s = tt R describes the junction with the braid (up to first order corrections in 
arc-length, as discussed below). 

6.1 Linearized Kirchhoff equations near a circular configuration 

For the loop solution, it is convenient to use the following cylindrical basis in the (y, z) plane: 



lerie) 



COS f e„ — sm f e» 



(45) 



1 eg (0) = sin Oey ~ cos 6 
These vectors defines an orthonormal frame (e,., e^, e^;) for any value of 9. With the choice 

m = |. 

this basis is adapted to the zero thickness solution in the sense that eg{9{s)) = t^(s). In the 
absence of ambiguity, the dependence of on s is not always written explicitly In the rest of this 
Section. The following derivation rules apply: 



de.r 
d^ 



eejO) 
R 



deg 

d^ 



erjO) 
R 



(46) 



Vectors decomposed in the basis (e^, eg, e^.) are denoted with square brackets. 

We introduce the first order perturbation of the tangent of the loop in the moving frame 
using two functions ii and v: 



ti(s) = t^(s)+e 



u{s) 


v{s) 



eg{s) + e (ti(s) er{s) + v{s) e^) . 



(47) 



By the inextensibility constraint (|2]), the perturbation to ti(s) along eg vanishes. The functions 
u and V are related to the functions xl, ijL and zl introduced in Eq. (32 1: 



= v{s), y'^{s) ^ ~cos{s/R)u{s), z'^{s) = - sm{s / R) u{s) . 
These equations will be used later to compute to x^, and z^. 



(48) 
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To use the constitutive relation (21 ), we first need to compute the derivative of the perturbed 



tangent given by Eq. (47 1. Using the derivatives of the cyUndrical vectors in Eq. (46 1, we find 



R 



u'{s)_ 
u{s) /R 



(49) 



Plugging this expression into the constitutive equation, we obtain 
mL(s) 



" 




' v'{s) + Uu{s) ' 




u 


+ € 


-v{s)/R 




l/R 


{r,e,x) 


-u'{s) + Uv{s) 





(50) 



an expression which is vahd up to first order in e. 

Like the tails, the loop is free of contact. As a result, the contact force vanishes, p = 0, 
and the internal force n(s) takes on a constant value over the whole loop. At dominant order. 



this value has to match that given in Eq. (28dl for the zero thickness solution. In addition, its 



linear correction in e has to be consistent with the symmetry conditions (24 1. This shows that 
the internal force in the loop is of the form 



R 



(51) 



where a and P are two constants to be determined. 



Combining Eqs. ( 50 ) and (511, we find that the equilibrium of moments ( 22c I can be expressed 



as a set of linear equations for the loop perturbation (m, v): 

R 

-u"{s) + Uv'{s) = ~(3 i 



R 



(52a) 
(52b) 



To integrate this differential system of total order four, we need to find four initial conditions. 



To this end, we proceed as in Eq. (51) and write down the centerline position, tangent and 



internal moment at the bottom of the loop which are compatible both with the zero radius 



solution, see Eqs. (28 1, and with the symmetry conditions (13 1: 

fL(0) 
tL(0) 



-2i?e, 



mi(0) 



1 

= 
R 



t/e^ + e (763; + (Se^) 



(53a) 
(53b) 

(53c) 



The constants p, (j), 7 and 5 introduced here will be determined later: p represents an infinitesimal 
motion of the bottom of the loop along the axis y of symmetry; 4) represents an infinitesimal 



rotation of the bottom of the loop about the axis y. Writing the two invariants, given in Eqs. (25) 
at s = we can eliminate to two other constants: 



J = (3 R and S - 



R 



(54) 



There remain four internal parameters for the loop, namely a and /3 introduced in Eq. (51 1, and 
p and (p in Eq. (53). Using a notation similar to that for the tails, we collect these unknown 
parameters into a vector "^l- 

*i - («, 13, p, 0) (55) 
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The initial condition for the set of differential equations ( 52 ) can now be written as a function 
of the loop parameters. They read 



u(0) = 


1 

e 


(ti(O)-t^(O)) -e. 


= 0, 


^'(0) = 


1 

e 


(mi(O)-mO(O)) •e, + [7«(0) 


= -I3R-U 


v{0)^ 




tL(0)-t"(0)) 


= -0, 


v'{0) = 


1 

e 


(mi(0)-m°(0)) •ey-C77i(0) 


= 



(56a) 
(56b) 
(56c) 
(56d) 



Equations ( 52 1 with initial conditions ( 56 ) are linear differential equations with constant 



coefficients. Their solution reads 

u{s) = a 5. i^Awf^ sin f 4 Kl{U)] - 4 ) (57a) 



Kim \Kl{U) \R ' 'J R 



R . (s 

sm 



|Xi(f/)^ {iiR+<i>U) 



-2 



«(.) ^ f c„ ( J) + „ =1=^ (cc, (I - l) (57b) 



where we have introduced the auxiliary function Kl{U) = {1 + R^ u'^Y^'^ . Integrating Eq. (48 1, 



one can find an explicit expression for the functions xl{s), VlIs) and zl(s) (the calculation is 
not difficult but the final expressions are long and the result is not given here). The constants 



of integration are provided by Eq. (53a) and are xl(0) = 0, ^^(0) = p and zl{0) = 0. 



6.2 Asymptotic expansion near junction with braid 

To match this solution with the braid, we shall need an asymptotic expansion of this solution 
near the top of the loop. In principle, this step is not difficult as it involves computing series 
expansion of the explicit solution just derived; in practice, the calculation is too tedious to be 
tractable by hand and was carried out with the help of a symbolic calculation language. 

It is convenient to describe the asymptotic shape of the top of the loop using a Cartesian 
equation. To do so, we eliminate the variable s in favor of z and expand the previous solution 
in series when s is close to tt R. Let us consider a current point on the centerline near the top of 
the loop with arc-length coordinate s — tt R + rj, where ij is a small quantity. We shall make a 
fundamental assumption, justified at the end, namely that rj is at most of order ^/e^. 



T]\ ^ e'^'\ (58) 



,1/2 
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To prepare the change of variable, we work out the relation between z and s 
z = zi(s) = z0(s) + ezi(s) + 0(e2), 

= + ^) + e 5l (tt i? + 7?) + ©(e^) , 

= -R sin (^^+ +ezL(7ri? + 77) +0(6^) 

= ^+ 0(fj3) + e (zi(7ri?) + 772^(71 i?) + 0(f )) + 0(e2) 
= 77 + ezL(7ri?) + 0(6^/2) 



In the last line we have used z'j\ti R) — 0, see Eq. (48 ). We have also collected all the O terms into 
a dominant contribution, of order 6"^/^ or smaller, using Eq. (58 1. Elimination of the arc-length 
variables s and 77 is then possible using the equality 



77 = s - TT i?^ = z - e 5l (tt i?) + 0(e3/2) 



(59) 



The term zi^{'k R) is known explicitly from the last section. Setting z = in Eq. (59 1 above 
yields the arc-length so of the point on the rod closest to the origin O: 



So = TT i? — e (tt i?) , 

where the quantity z^ (tt R) in the right-hand side will be given at the end of Section [s] 
We expand xi^ and yj^ similarly to z^: 

xl{-kR + v) = e {xLiTrR)+rjx'L{TrR))+0{e'^) 

(vr i? + Tj) = - ;J + e (yi (^ i?) + 77 y^TT i?)) + O(e') . 



(60) 



(61) 



The first term in the right-hand side of the second equation, —rf/{2 R), arises from curvature of 
the loop in the zero thickness solution. This term has to be retained as it is of the same order of 
magnitude as the other terms when 77 is of order ^/e. 

We can now use Eq. (59 1 to eliminate the arc-length 77. This leads to a Cartesian equation 
of the top of the loop in the form: 

XL{z)^eXL + ezX'^ + Oie^), 



2R 



eYL + ezY[ + 0{e^), 



Because of Eq. (58 1, this expansion is valid when z is of order ^/e or smaller: 



1/2 



(62a) 
(62b) 

(62c) 



The coefficients {Xl, X'j^,Yl,Y[) of the asymptotic expansion are found by identification with 
Eq. (leTl): 



XL = XL{nR) 
X'l^x'l{t:R) 



YL = yL{T^R) 
Y[^y^{7TR) + 



zl{^R) 
R 



The quantities in the right-hand side are all known explicitly from the analysis of the loop given 
in Section 6.1 Being solutions of a set of linearized differential equations, they all depend linearly 
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Figure 5: Braid geometry. We call a the strand having positive arc length s (s « ttR in the 
center of the braid), and b the strand with negative arc length (s « —ttR in the center). 



on the loop parameters '4'l = {a, /3, p, (f>), as revealed by Eqs. (57 1. The linear mapping giving 
the expansion coefficients as a function of the loop parameters reads 



Yl 

V n ) 



Mz,(;7) • *L, where Ml({7) 



V 



(Kl/W 
1+2B?U'^+Kl 
UKl 



UKl 
_(l + n) 
U_j_R 

~iFJr~ 



_u _ 

RU 



\ 



(63) 



To keep the notations compact, we have noted Kl = Kl{U) = (1 + i? U Y^"^ and introduced 
the shorthand notations = 'K^iU) =_cos{TrK l(U)) and = kI{U) = sm{Tr Kl(U))- This 
explicit expression for the matrix Ml(C/) comes from the an alyti cal solution for the loop given 
in the previous section. Recall also that R = 1/V2 by Eq. pol. The matrix Ml{U) defined 
above PI plays a role similar to M.t{U) for the tails: it captures the elastic response of the loop 
to the perturbation induced by the presence of the braid. 



7 Braid solution 

The solutions in the outer domain (loop and tails) have been derived in the two previous sections. 
We now proceed to solving the internal region (braid), which is more difficult as it involves self- 



contact. A key remark is that the scaling relations expressed in Eq. (33 1 imply that the tangent 
deflects from the z axis by a small angle, of order h/I e <^ 1. As a result, the approximation 
of small displacements hold and the Kirchhoff equations can be linearized; by linearity, the braid 
problem can then be decomposed into an average problem without contact, and a difference 
problem where contact takes place with a fixed, virtual cylinder. As earlier, a subscript B 
denotes quantities associated with the braid. The two strands composing the braid are labeled 
with superscripts a and b, as shown in Fig. |5] 

7.1 Centerlines 

We have introduced a rescaled axial displacement consistent with the scalings for the braid: 

<j=- = ^ (64) 



Note that the matrix M.t{U) has a smooth hmit for C/ — > 0. Even though there are some powers of U in the 



denominators, the following expressions are smooth near U = 0: and =2- - J 
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The leading order term of the expansion in the braid was given in equation ( 33 1 . For the first 
strand, labeled a, it reads 



The centerline of the other strand is defined by a similar formula with a replaced by h (note that 
the constant tb, which represents a global translation of the braid, is common to both strands 
and so has no index a or b) . Our unknowns for the braid problem are the translation Tg and the 
four functions x^, ijg, x^g and y^, defined in terms of stretched coordinates. 

Recall that the strands a and h are mapped onto each other by a symmetry of angle tt about 
the z axis. By our choice of axes, tr = is the center of the braid: the symmetry is expressed by 
the following relations, 

^s(^) - -S^b{-^) (66a) 
Vsi^) = +rB{-^), (66b) 

which implies that there are only two independent functions to be determined, say x°' and 

It is useful to introduce an auxiliary quantity, the velocity d^io") at which the centerline is 
swept out in this parameterization — we do not use arc-length parameterization here: 

c'^{a) = rB{a)\ = e + 0{e^). 

The unit tangent is then defined by 

^b{^) = ^IS = + ^ ^^b\^) + VBi^) e.) + • • • (67) 

C°((T) 

and a similar equation for the other strand. Note that primes applied to functions such as xg, 
ys or denote derivatives with respect to their argument, 7f here. 



x%{a) 
e a 



(65) 



7.2 Contact 

For any pair of points in contact in the braid, let ct" and be the rescaled coordinates of the 
point on braid a and on braid b, respectively. The contact condition writes 

\r%{a-)-f''g{a')\^'^ = V2e\ (68) 

We square both sides of the equation and use the the centerline parameterization given in equa- 
tion (65): 

[{x%{an - xU^"))' + {yU^n - y'si^'m + {a^ - a")' = {V2e'r 

In this equation, there is only one term of order and no lower order term; this term has to 
cancel, which implies: 

ct" — a'' (for points in contact) (69) 

At next order, we obtain 

{x%{a) - x''g{a)f + {y%{a) - i/gia))^ = 2 (for points in contact) (70) 
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By equation (69), contact occurs only between points lying in the same plane perpendicular 



to the z axis at the leading order in e. Let us define the locus of the contact in physical space as 

£) {ct such that |r^(a) - r^(CT) | V2 e^} 

This set £) is composed of the rescaled z coordinate, called (f , of the points in contact, unlike the 
original set £ which describes contact points based on pairs of arc-lengths. In the limit e <C 1, 
the set S) provides a description of contact much simpler than the generic one, based on £. In 



Eq. ( 88 ) , we shall compute the set £) explicitly, and show that it has non-trivial topology (it has 
'holes' in it). For reference, we mention that the initial set £ can be reconstructed by 

U {(si,S2),(s2,Si)} where y {((so + e^)L*, (-so + e'^)i'')} 



{si,s2}e£' 



and So was defined in Eq. (60 1 as the arc- length of the point on strand a closest to origin, that 
is such that z = 0. 

We s hall now express the contact pressure p. First note that the action-reaction principle in 
Eq. (Ill can be rewritten as p°'{a) = p^((t); we can therefore omit the superscript and note p{a) 



the scalar contact pressure associated with contact occurring at coordinate z = ea. According 
to equation (|9|, the (vector) contact pressure is the scalar pressure p{a) times the unit vector 
joining the barycenters of two cross-sections that are in contact; in rescaled form this reads, see 



Eq. (201 



p''(a)=p(a)^"«^^) 



-b -r—\ I x%(a) — X^^((t) 





This quantity appears to be orthogonal to the z axis at this order, as expected. 



(71) 



7.3 Equations of equilibrium at leading order 



Combining the constitutive relation (21 1 with the formula (67 1 for the tangent, we obtain the 
internal moment as: 



c°(ct) 




(72) 



Note the normalization factor l/c° in the above expression for the normal curvature vector, 
dt/ds = {di/d(j)/c"-, which is required since parameterization does not use arc-length. 

For any vectors a and u such that u has unit length (u^ = 1), the following identity holds 
a = (a • u) u — u X (u x a). With a ~ ng((f) and u — t°g{a), it can be used to compute the 
internal force: 

n%ia) ^ {tlia)-n%{a))tl{a)^tlia) x (t^(^) x n%{a)). 
The factor (^^(ct) • H^ct)) can be expressed using the second invariant I2 = U /2 + 1 given in 



Eq. ( 25b I , while the balance of moments ( 22c I allows one to rewrite the vector in the last term 
as (t^(CT) X n^(CT)) = —'^^Sj=y-, this right-hand side being given itself by Eq. ( [72| . This yields 
the following expression for the internal force in the braid 




m. 



—a , , —a , , 

ti3(CT) + tB(cr) X 



m|'(a) 

c'^(ct) 



(73) 
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The first term in the right-hand side is of order (it is bounded for small e) and is dominated 
by the second term, of order 1/e because of the denominator c°'{a) — £+•••. This yields the 
leading order term for the internal force: 

where the ellipsis stands for negligible terms that are bounded for small e. The internal force 
in the other braid is given by a similar formula. 



The balance of forces (22dl writes n%' (a) / c"" (a) +p(ct) = in the current parameterization. 



Using Eq. (71 1 for the contact force, this yields 



and similar equations for the y direction and for the other strand. This equation shows that the 
rescaled contact pressure p has to be of order in order to balance bending. Therefore, we 
define the final rescaling for the contact pressure by: 

g2 ^1/2 

p{a) = e^p{a) = ^^^^ p{a). 

By the previous argument, this p(jf) has a finite limit for e — > 0. The equations of equilibrium 
then take the form 

x%""{a) = ^p{a) {i%ia) ~ 4(a)) (75a) 

y%""{a) = ^p{<j) (y^(a) - y|(a)) . (75b) 

This is a set of fourth-order differential equations for the deflection, which arc coupled through 
contact. The contact pressure p{a) is the Lagrange multiplier associated with the non-penetration 
condition, and is not known in advance. The rest of Section[7]is devoted to solving these equations 
with appropriate boundary conditions. 



Equations (75 1 stand for an elastic rod in the small deflection approximation, subjected to the 
normal distributed force given by the right-hand side. Note that the twist loading parameter U 
does not appear in the equations for the braid at the leading order, as bending effects dominate 
twist. A nice consequence is that the braid problem is universal: unlike the tail and loop problems 
studied earlier, the formulation of the rescaled braid problem involves no parameter (except for 
the knot type which is a discrete parameter). 

7.4 Decomposition into average and difference problems 



Taking advantage of the linearity, we can combine Eqs. (75 1 and the two similar equations for Xg 
and into a difference and an average problem. The average variables (/, g) and the difference 
variables (u, v) are defined as follows: 

f{a) = {x%{a) + xU^)) u{w) = ^ {x%{a) - x%{w)) (76a) 

g{<f) = (y^(a) -I- y^(a)) v{a) = {y%{w) - y^(a)) (76b) 
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Summing Eq. (75a I and the similar equation for strand b, namely 



b nil 



we 



find /"" — 0. By the same argument, g" 
average problem: 



0. The unknown contact force disappears from the 



/""(a)=0 
5""(a)=0. 



(77a) 
(77b) 



In the next section, we derive the asymptotic conditions associated with these equations, which 
are then solved in Section FTHl 

For the difference problem we obtain 

u{a)"" = {^f2v{p))u{a) (78a) 
v{a)"" ^{^v{o))v{a) (78b) 

Although the contact force p(ct) is still present in the difference problem, the contact condi- 
tion ( 70 ) takes a very simple form when formulated as a function of the difference variables: 
V? + ■ir = 1. As a result the non-penetration condition is expressed by the inequality 

u'{a)^'V^{a)>\, (79) 

and the problem is much easier to solve. The average and difference variables are subjected to 
the following parity conditions, deriving from equation (66 1: 



/(-a) = -/(a) 
5(-^) = 9^) 



u(— cr) = u((t) 



(80a) 
(80b) 



7.5 Asymptotic conditions 



We have derived in the previous section the differential equations for the braid. These equations 
make use of a stretched variable a. In the present section we derive the asymptotic conditions 
the solutions must satisfy for large values of the stretched variable a. As usual in matched 
asymptotic analysis, the asymptotic conditions for the inner problem are required for the inner 
solution to match the outer solutions in the region of overlap (intermediate region), where both 
the inner and outer solutions are valid — see Section [8] for a detailed discussion of this matching 
procedure. Given our conventions, summarized in Fig. [5] the strand a of the braid connects to 
the loop for ct — > — cx), and to the tail for tf — > -|-oo. 

Let us start with the condition for matching the internal moment. At the top of the loop ~ tt 

and so ~ e^; the internal moment, given by Eq. (50), then reads mL(7ri?) = e^^/i^-l-f/e^-t- 

Comparison with the braid moment, given by Eq. (72|, yields the asympt otic condition x°^' ~> 
and iIq" ' - - 



1 



oo. Using the value of i? = l/-\/2 given by Eq. (30), we write 







for a 



Vb 



—a/2 for (T ^ — c 



The asymptotic condition f or a — > -|-oo is obtained by matching the moment at the origin of the 
tail, m.T = UBz, with Eq. (72). It yields: 



Vb 



for (T 
for a 



'OO, 

-oo. 



26 



A similar argument gives the matching condition for the internal force. The force is bounded 



for small e in the two outer regions, see Eqs. (35) and (51 1. In contrast the force in the braid 

i. Tt 

-oo): 



diverges for small e, as shown by the term of order (1/e) in Eq. (74 1. This term must therefore 
vanish when strand a reaches the loop {a 



-oo) or the tail (cr 



x%"' for a ±oo, 
yy for a ±oo. 



Using the parity conditions (66 1, we obtain the same equations for the other strand. The asymp- 



totic conditions arc then expressed in terms of the average variables: 



n±^) - 

/"'(±oo) ^ 



and of the difference variables: 



u"(±oo) ^ 
w"'(±oo) ^ 



<7"(±oo) ^ -1 
g"'(±oo) ^ 



v"{±oo) Tl 
v"'{±oo) ^ 0, 



where we used a shorthand notation meaning that v" goes to —1 for a 
a — > — oo. 



(81a) 
(81b) 



(82a) 
(82b) 

-oo, and to +1 for 



7.6 Solution of average problem 

The average problem being insensitive to contact forces, its solution is straightforward. The 



general solution of Eqs. (77 1 yields for /((t) and g{a) polynomials of order 3. To satisfy the 
parities, see Eq. (80 1, we write: 



g{a) = Co + C2CT^ 

where cq, Ci, C2 and C3 are real constants. Two of these four constants are set by the asymptotic 



conditions (81 1 and we have: 



/(ct) = Ci cr, 

9{^) = Co - 



(83a) 
(83b) 



The two remaining constants co and ci will be found by solving the matching problem, see 
Section [1 

7.7 Solution of difference problem 

We proceed to solve the equations for the difference problem (|78|, subjected to the asymptotic 



conditions (82 1. In the right-hand sides of Eqs. (78 1, the unknown contact pressure p{a) has to 
be determined, in a way that is consistent with the contact set D: pijr) can be nonzero for those 
(f that are elements of D only. By Eq. (79) the contact set J) depends on the difference variables 
only: 

D = {ct such that u^(ct) -I- v'^{a) l}. (84) 
This set will be found as an outcome of the solution of the difference problem. 
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7.7.1 Variational formulation 



Solution of the difference problem is greatly eased by pointing out the simple variational structure 
underlying the equations. We shall now show that solutions of the difference problems are 
minimizers of the following energy: 

!M^^)+^da + .W + z;'(-W^). (85) 



J-w ^ 

Minimization is done with respect to the functions u{(j) and v{7j) which are defined over the 



interval [— VF, W\ , are twice differentiable, and are subjected to the non-penetration condition ( 79 ) 
— here, W is large but fixed number. We shall also include an additional constraint, related 
to the knot type: the parametric curve (u((f), ^(fj)) has to make a prescribed number of turns 
around the origin; this topological constraint is discrete and so does not affect the Euler-Lagrange 
equations. 



Before we show that solutions of the braid equations (78 1 subject to asymptotic conditions 



and constraints are minimizers of the energy (85 1 for large enough real numbers W , we shall 



first give a physical interpretation of this energy. To this end, we define a virtual rod, called 
the difference rod, by the following Cartesian equation: \x = u(j7),y = v(jf),z = (i/e}. Being 
based on the difference variables u and v, this difference rod winds around the z axis exactly 
in the same way as the strand b winds around the strand o in the original problem. Note 
that this difference rods extends from z — ~W/e to z = +W/e: because of the factor 1/e in 
the definition the rod deviates only slightly from the z axis and its arc length is approximately 
z « fj/e. Using the small deflection approximation for this difference rod, one can easily compute 
the unsigned curvature of its centerline, n = {u" + v" )^/^. By definition, the difference rod 
has zero twist, and we define its bending modulus to beSdiff = 1- Then its elastic energy ([4| 

reads J^^!^ ^ dz — - — — — da. Up to the factor e'^ , which is irrelevant for the 



minimization problem, this is exactly the first term in our energy (85 1. The two remaining 
terms, v'{±.W), can be interpreted as the work of the bending moments on the endpoints: the 
unit tangent to the rod is {eu',ev', 1} and so a moment (e^e^;) applied on either end of the rod 
is associated with the potential energy e^v'{ztW). The quantity is then factored out of the 



total energy. The constraint ( 79 1 is interpreted by the fact that the difference rod winds around 



a virtual cylinder whose axis is the z axis, with unit radius. To sum up, we have identified 



the energy (851 as that of a virtual, twistless, naturally straight rod winding around a fixed 
cylinder, of unit radius and axis e^, and subjected to bending moments at its endpoints. We 
have transformed the self-contact problem into a contact problem with a fixed, external body, 
and this an important simplification. 

We shall now establish the equivalence of the constrained minimization problem and the 



original braid equations ( 77 1 , by working out the Euler-Lagrange equations for the minimization 



problem. First, let us rewrite the constraint (79) as Q > 0, where 



"^^"^^ — 71 71- 

Constrained minimization problems are classically solved by introducing Lagrange multipliers, 
here the function 7r((7), and enforcing stationarity of the augmented energy. 



r+W 

SE- n{a)SQ{a)da ^0. 

J-w 
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Using the explicit expressions of E and Q given above and integrating by parts, this yields: 



r + W 

-w 



u du + {v"" - V2 TT{a) vj Sv 



da 

+w 



-w 

+ {v"{W) + 1) dv'{W) + {-v"i-W) + 1) Sv'{-W) = 0. (86) 

Here, square brackets with subscript and superscript denote boundary terms coming from the 
integration by parts, [f]^ = f{b) — f{a). The quantity in the left-hand side has to be zero for 
arbitrary variations Su(a) and 6v(a). Therefore, the factors in front of Su(w) and Sv{a) in the 
integral have to vanish: after identification of the Lagrange multiplier ^{a) with the rescaled 



contact pressure p, one recovers the equations ( 78 ) of the difference problem. The remaining 



boundary terms in the variation above yield the following boundary conditions: 

u"{±W) = v"i±W) = Tl (87a) 

u"'{±W) = v"\±W) = 0. (87b) 

For large [^VT, we recover the boundary conditions (82) derived earlier for the difference problem. 
This establishes the equivalence of the two formulations. 

7.7.2 Numerical solution of the universal braid problem 

We have just reformulated the difference problem as a constrained minimization problem. We 
now take advantage of this variational formulation and present a numerical solution which is 
very easy to implement. The difference problem has been formulated without any parameter: 
for any given knot type, the solution of the braid problem is universal. In particular, note that 
the twist parameter U has been removed from the braid equations at dominant order: the braid 
is insensitive to the applied twist. These universal solutions are computed below, once for all, 
for the trefoil (3i) and cinquefoil (5i) topologies. 

We first implemented the minimization problem using the symbolic calculation language 
Mathematica which has built-in capabilities for non-linear constrained optimization. The imple- 
mentation is straightforward. The problem is first reformulated in polar variables {w,(j3), such 
that u — w cos (f) and v — w sin 0: the advantage is that the winding number about the z axis 
is readily available from the end value of (j). Values of the functions w and are sampled on a 
uniform mesh covering the positive axis, a g [0, W], and their values for negative a are recon- 



structed using the parity condition (80). Finite differences are used to evaluate the objective 
function ( |85[ ). The non-penetration condition ( |79[ ) is enforced by a constraint w > 1 written at 
every point of the mesh. In addition, we use a series of non-physical constraints: (i) we require 
that (n — i) TT < (f>{W) < {n + 1) tt, where n = 1 for a trefoil knot and n = 2 for a cinquefoil 
knot; (ii) we require that \4>{ai) — 0((Ti+i)| < ^ for any pair of neighboring mesh points ai 
and (f,j+i. Constraint (i) is used to direct convergence towards the solution having the required 
winding number, as the difference rod has to make one and a half turn around the z axis in the 
trefoil case and two and a half turns in the cinquefoil case — note that the winding number is 
given by {(t){W) — 4'{—W))/{2t:) — 2 0(VK)/(27r). Constraint (ii) warrants that cj), defined mod- 
ulo 27r, varies smoothly along the rod which is required for the end value (p{W) to express the 



^Convergence of our variational problem for large W is extremely simple: as we shall show, the minimizer 
becomes independent of W when W is larger than a fixed number, which can be interpreted as the coordinate of 
the last point of contact with the cylinder. 
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knot type 




(Tp 




A 


3i 


0.348 


2.681 


1.823 


0.022 


5i 


4.504 


6.814 


5.962 


0.021 



Table 1: Numerical values of the contact-set parameters for 3i and 5i knots. 



total number of turns. We carefully checked that the non-physical constraints (i) and (ii) are 
non-active when the minimization procedure exits, i. e. all inequalities are strict: their role is 
simply to guide convergence towards a physically relevant solution. 

We found that the minimization always converges to the same type of solution for both knot 
types, n = 1 or n = 2. We used a typical mesh size of Act ~ 0.1 and interval width W ^ 9 — 
we observed that the numerical solution does not vary with W when W becomes larger than 4, 
something that we shall explain soon. In Fig. [6] the difference rod is visualized for the trefoil 
topology. By inspecting where the constraint w > 1 is active in the numerical minimizers, we 
can determine which mesh points belong to the contact set D. When the mesh is not exceedingly 
coarse Act ^ 1.5, and for both knot types, we found an interesting contact topology: the contact 
set is composed of an interval centered around the origin and two symmetric isolated points (each 
corresponding to a single mesh point). Starting at ct = 0, the difference rod is in continuous 
contact with the cylinder, then lifts off from the cylinder and eventually touches it again at an 
isolated point. This contact set is shown in Fig. |6] Note that this topology remains the same 
when the mesh size is decreased. This leads to the following topology for the contact set D of 
the difference problem: 

D = {-CTp} U [-CTe, CTe] U {Wp} whcrC < CTe < CTp. (88) 

Here is half the width of the central region with continuous contact and CTp is the rescaled 
coordinate of the isolated point of contact. We stress that this topology arises from the numerical 
minimization without any a priori assumption on our part. In a problem where contact occurs 



along straight line in space, Coleman and Swigon (20001 have assumed a topology of this form 



and checked that it was consistent. 

By our definition of the difference variables, the distance w of difference rod to the z axis is 
also the rescaled distance between the centerlines of the strands a and b in the original problem: 
contact of the difference rod with the virtual cylinder {w = 1) means that the physical strands 
a and b are in contact. Like the virtual bodies, the two physical strands experience continuous 
contact in a central region; on both sides of this central region, they separate by a small but 
finite distance, contact again at a point, and finally separate for good. The maximal reopening 
A is given by the extremum of the function (w — 1) in the interval [CTe,CTp], see Fig. |6j the 
corresponding value of ct is called Wg. These numerical values are given in Table [T] Values of A 
are very close for trefoil and cinquefoil knots, A w 0.021; in physical units, this corresponds to an 
inter-strand reopening of (0.043 h), that is 43 fim for a rod of radius h ^ I mm. The experiments 
reported in Section [9.2| confirm the presence of these openings. 

In order to confirm our hypothesis on the topology of the contact set, we implemented an 
independent numerical solution for the difference problem, assuming a topology of the form 



This independent solution relies on non-linear shooting: in contrast to the energy minimization 
scheme, it involves a numerical integration of the equations of equilibrium; this new approach 
is much more accurate but requires the contact topology to be known. Numerical integration 
is carried out on each interval ct G [0;CTe], [CTe;CTp] and [CTp;cx3] in turn. In the first interval the 
difference rod lies on the surface of the virtual cylinder with unit radius, and we use the polar 
variables {w,(j)), introduced earlier, with w = 1. The polar variable 4>{a) satisfies the differential 
equation 0"" = 6{(f)')^ 4>" . Four initial conditions are required to integrate this equation, two of 
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Figure 6: Numerical solution of the difference problem of the braid for the trefoil topology 
(n = 1). The difference rod, shown in red, describes position of strand b with respect to strand 
a — compare with Fig. [5] It is held by bending moments applied at its endpoints, and enlaces 
an impenetrable cylinder of unit radius drawn around the z axis. The problem has no parameter 
and the solution depends on the knot type only, (a) 3D view, (b) projection onto the plane 
(u, v) perpendicular to the cylinder axis, (c) distance w = ^/u'^ + v"^ of difference rod to cylinder 
axis: w = 1 when there is contact, and w > 1 otherwise. The contact set is denoted by shaded 
regions (in blue) along the solution. Note the interval of contact around the center of symmetry 
(— (Tg < (7 < (7e), flanked by two isolated points of contact {a — iijp). Close examination of (b) 
reveals reopening in the intermediate regions (ffg < \a\ < Wp). 
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which are fixed by the symmetry condition (80), 0(0) = and 4>"{0) = 0; the other two, 0'(O) and 
4>"'{0), are unknowns of the shooting procedure. All the other quantities can be reconstructed 
from </i((t). At Ife a jump Pg in the internal force is introduced. It represents a Dirac contribution 
to the contact pressure The values of 4> and its derivatives at the end of the first interval 
are combined with P. to evaluate the initial conditions for the second interval. In the second 



interval, there is no contact and Eqs. ( 78 1 are integrated with p{a) = 0. At Wp the rod touches the 
cylinder and there is another discontinuity Pp in the internal force. In the last interval [ap] +oo] 
there is no contact and the internal force is again constant. By the asymptotic conditions (82) 
this constant force has to be zero. This implies in turn that the internal moment is constant. In 



view of this the four asymptotic conditions (82 1, which concern the internal force and moment, 
have to be be satisfied over the entire third interval. Overall, the shooting scheme involves six 
unknowns {0'(O), 0"'(O), (fe, (fp, Pg, Pj,} which must satisfy six equations, namely two geometric 
contact conditions at ap and four conditions coming from Eq. ( 82 1 . For the trefoil knot the 
non-linear shooting procedure converges to {0'(O) = 0.769, 0"'(O) — 0.033, = 0.348, CTj, = 
2.681, Pg — 0.170, Pp — 0.442}. The contact pressure can be reconstructed in the first interval 
as Tr{a) = p{a) = - 30"^ - 4(/)' (j)"')/V2. It is plotted over the full contact set J) in Fig. [tJz: 
it is everywhere positive and this validates our assumption on the topology (in Appendix |A] we 
test different contact topologies and show that they lead to negative pressure and/or residual 
penetration) . 



The internal force is given by Eq. ( 74 1 for strand a and a similar equation holds for b. Noticing 



that the third derivatives of the average solution given by Eq. (83 1 vanish, we find the nonzero 
components of the internal force in each strand: 



(n^). 



U"'((7) 



v"'(a) 



(89) 



The rcscaled internal force is plotted in Figure |7|). Note the discontinuities of the internal force 
at the boundaries of the contact set, where Dirac pressure forces are present. 



7.7.3 Polynomial expression beyond last contact point 



In Fig. [7] the rescaled internal force appears to be zero beyond the isolated contact point, that 
is for |(t| > (Tp. From Eq. (89), the third derivatives of the functions uijf) and v{a) vanish in 
this region. As a result, both u and v are polynomials functions of a of order at most three. 



In addition their cubic term has to be zero for the asymptotic conditions (82a) to be satisfied. 
Therefore both u(ct) and v(w) are second order polynomials for |ct| > ap. The quadratic term is 



fixed by the other asymptotic conditions (82a I: it is zero for u{a), which is therefore an affine 
function, and it is — | for w((t). Consequently, for (f > (fp, u is of the form u{(f) = A„(f-|- qn for 
some real constants A„ and and v is of the form vijf) = — ^+n„ (t+(7^ for some constants n„ 



and q'„. The expressions for ct < — (fp are found using the parity conditions (80). The following 
condensed notation summarizes both cases a > Up 
as ±(7 > (fp): 



and <j < —Gp (which are denoted generically 



u{(t) 
v{a) 



±Ar, cr - 



(90a) 
(90b) 



In this condensed notation, one should use the upper sign on the positive side, for +cr > cjp, that 
is replace ± with (+) and =F with (— ), and the lower sign on the negative side, for —a > Wp. 

^"Dirac contributions to the contact pressure ap pear generically at the boun dar y of the contact set in contact 
problems for clastic rods, as shown for instance by Coleman and Swigon 1 2000 I, or Audoly and Pomeau (2008i. 
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Figure 7: Forces in braid for the trefoil geometry, same solution as in Fig. [6j (a) rescaled contact 
pressure; (b) rescaled internal force. The pressure is everywhere non-negative and this validates 
the assumption on the contact topology. The internal force is proportional to m'" and w'" by 



Eq. ( 89 ) . The localized contact forces Pe and Pp are represented by columns in (a) , and manifest 



themselves as jumps in (b) . 
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knot type 


An 


n„ 


3i in=l) 


-0.87759 


2.089 


5i (71 = 2) 


-0.87738 


6.223 



Table 2: Numerical values of braid constants n„ and A., 



The coefficients A„, n„, g„ and are available from the numerical solution of Section 7.7.2 
The values of A„ and n„ are given in Tabled As we shall see later, the values of g„ and are 
irrelevant at dominant order and are not given here. 



7.8 Asymptotic expansions at braid-tail and braid-loop junctions 

In order to match this inner solution with the outer solutions computed earlier, we shall need its 



expansion far away from the braid, that is for large values of cr. In Eq. (76 1 the inner solution 



is decomposed into an average and a difference solution. The average solution is a polynomial 



given by Eq. (83). The difference solution is polynomial as well for large enough values of a, see 



Eq. ( 90 ) . It is straightforward to combine these polynomials to obtain the expansion of the braid 
solution for a — > ±00: 



V2 



Here, we use the same condensed notations as in Eqs. (90 1, whereby the compound signs ± and 
=F must be replaced by the upper symbol when a —^ +00, and by the lower one when a —^ —00. 
Note that we have replaced the stretched variables x% and a with the barred variables x% and 
z using Eq. ( 65 1 , as the matching has ultimately to be done using a common set of variables for 



the inner and outer solutions. 

A similar calculation for yg yields: 



V2 



(g(cr) - v{a)) 



(^(-l±l)^-n„a+(coT9;) 



z 

71 



-1 ± 1 



e \TB- z-j^]+0{e'). (91b) 



In this equation, the coefficient tb represents an infinitesimal translation of the braid along the 



y axis. The term proportional of ci in Eq. (91a I is very similar: it represents an infinitesimal 



rotation of the braid about the y axis. We rename it ws, 

lob 



V2- 



(92) 



The two other coefficients in Eqs. (91 1 have been computed in Table [2j We call internal param- 
eters of the braid the two remaining free parameters in the above expansions: 



*B = (ts, lob)- 



(93) 
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In the next section, we shall show how these parameters '5'b can be found as a function of the 
applied loading and knot type, together with the loop and tail parameters 'S'l and \I't- 



We can rewrite the expansions (91a I and (91b I in the form 



eX 



-1 ± 1 



z 

2^ 



(zX 



ezY, 



B 



(94a) 
(94b) 



where i? is a shorthand for 1 / \/2 by Eq. ( 30 ) . Note that the factor in parenthesis in Eq. ( 94b I is 
equal to ( — 1 — l)/2 = —1 on the negative side, and to (—1 + l)/2 = on the positive side. As a 
result, the quadratic term in y%{z) is equal to — z^/ (2 R) for large negative a, which is consistent 



expansion (43b I for the tail 



with expansion (62b I for the loop, and is absent for large positive cr, which is consistent with 



The coefficients of the polynomial expansions just written are given by identification with 



Eqs. (911: 



xf 



Y, 



= M± • *B + V± (n) , where M± = 



r 







1 


1 







0/ 



and V^(n) 







(95) 



Equation ( 95 ) defines two constant matrices and , and two vectors [n) and [n) 
depending on the Icnot type n. These vectors are defined in terms of the braid constants found 
in Section ITTZl 

For the matching problem studied in the next Section, it is useful to give a precise description 



of the range of values of z where the expansions (94 1 hold, that is where the omitted terms 
denoted by ellipses are actually negligible. These expansions have been obtained by taking 
the limit \a\ — s- oo: they obviously require \a\ 3> 1, that is \z\ ^ e. However, this is not 
the only assumption. Recall that the braid has been studied based on the small displacement 
approximation, which assumes that the tangents or tg remain close to the vector — see 
for instance the tangent expansion (67 1. This assumption breaks down in the inside of the loop. 



for values of z of order 1: as shown by the quadratic term in Eq. (94b I or directly by the loop 



solution (28b I, the tangent deflects from the z axis by an angle {z/R) there. Therefore, the braid 
solution accurately describes the upper part of the loop, where it merges with the braid, but 
does not accurately describe the whole loop: it assumes |z| ^ 1. To summarize, the range of 



validity of the braid expansions ( 94 1 is 



e < \z 



< 1. 



(96) 



The linear relations in Eq. ( 95 ) yield the braid expansions in the regions where it connects with 
the loop (—1 ^ z ^ — e) and with the tail (e <C z <C 1). These relations depend on the internal 



parameters of the braid, "^b ^ (tb, ^b), and on the knot type n. These equations (94 1 and (95) 
capture all what we need to know about the inner solution (braid) to be able to solve the problem 
globally. 



^^As explained in Section |7.7.3[ the braid actually reaches its asymptotic behavior exactly as soon as the 
last point of contact is passed, |(t| > (fp. This is not important and we shall continue to write the less severe 
requirement \a\ 3> 1, which holds in general in boundary or inner layer analysis. 
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8 Matching 



So far, we have solved the equiHbrium equations in the tail, loop, and braid regions independently. 
In each region the solution depends on some parameters, collectively denoted 'U't, '^l and ^b, 
which have yet to be computed. Figure |4] illustrates the fact that these domains overlap in the 
so-called intermediate regions. There are two types of intermediate regions, one where the loop 
merges with the braid, and the other one where the tail merges with the braid. By writing down 
the matching condition of the various pieces of solutions obtained so far in these intermediate 
regions, we make sure that we have constructed a smooth, global solution of the original problem. 
We now derive these matching conditions and compute the remaining parameters \I't, l and 

8.1 Matching braid and tail 



In the end of our analysis of the tail regions, in Eq. (43 ), we have obtained the following expansion: 



xt{z) = tXx + ez X'rp 



e z Yr^ 



(97a) 
(97b) 



where the terms that have been dropped, of order and ez^, are negligible if z <C e^^^. In 



Eq. ( 94 ) , we have found a similar expansion based on the braid solution: 



x%iz) 
TbC^) 



eX+ 



-ezX+'- 
e-zY^'^ 



(98a) 
(98b) 



which is valid for e ^ |z| ^ 1. These two expansions have to be consistent in the region of 
overlap, defined by e ^ z <C e^^^, and this implies the equality of the coefficients. We obtain the 
following matching condition in the intermediate region between braid and tail: 













Yt 




n 


Vt) 







Using the reduced matrices and vectors of the tail and braid problems defined in Eqs. ( |44[ ) 
and ( 95 ) , this matching condition is rewritten as a linear system for the tail variables = (A, /z) 



and the braid variables b — {'''Bt'^b)'- 

Mt{U) • *t = M 



which reads: 



/ 1_ 

-a{U) 




-biU) 

1 

-a(U)J 



fo 0\ 

1 

1 
\0 0/ 



V+(n), 



frB 





V2 



(99) 



The functions a{U) and b(U) were defined in Eq. (39). 

This is a set of four linear equations for the four unknowns A, fi, tb and ub- As can be 
checked easily, the determinant of this linear system is a{U). For U ±2, a{U) ^ and this 
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system has a unique solution, which can be found exphcitly by ehmination: 



X{U, n) = 



tb{U, n) 



ujB{U,n) 



V2 



(100a) 
(100b) 



after using the detailed expressions for a{U) and h{U). 

We have just found the internal parameters of the tail and of the braid, and ^ b, as a 
function of the dimensionless parameters of the problem, the loading parameter U and the knot 
type n — recall that the braid constants n„ and A„ were given in Table [2] for trefoil (n ^ 1) and 
double knots {n — 2). 

8.2 Matching braid and loop 

A similar argument holds for the intermediate region between braid and loop. In Section [5. 2[ we 
found that the top of the loop is accurately described by the expansion 



z 



eYL + ezY[ 



(101a) 
(101b) 



This expansion is accurate up to terms of order e^, provided {—z) <^ e^^^, and of course z < 0. 
On the other hand, the braid solution has been expanded as 



x%{z) = eXg + ezXg' H 



(102a) 
(102b) 



in the domain defined by — 1 ^ ^ <C — e. The intermediate region between braid and loop is 
defined by — e^^^ <C z <C — e. There, the two expansions have to be compatible, which leads to 
the matching condition: 



(XA 

X'l 




_ Xb 


Yl 






WJ 



As earlier, we arrive at a linear system which can be written in terms of the reduced matrices 
and vectors of the loop and braid problems, defined earlier in Eqs. (63 1 and (95 1: 

Ml(U) ■^L = Mg -^B+ Vs(n). 

We arrive at a linear system for the loop parameters: 



P 







\/2 



tb 
n„ 

72 



(103) 
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where the quantities tb and ub in the right-hand side are given by Eq. (100b) and the 4x4 



matrix M.l(U) is defined in Eq. (63 1. The determinant of the matrix Ml(C/) can be computed 
exactly; it vanishes for 



detML(;7) = iff U = ±^2 {f - 1), for some integer j > 2. 



(104) 



For any other value of U , one can solve the linear system in Eq. ( 103 1 by elimination. This yields 
the following expressions for the internal parameters of the loop: 

a(U,n) = 

TT 

17 n„ uV2ii„ u 



PiU,n) = U\ -2A„ + 



V4^ 7r(2 + L/') y~=^tan(^if,.(C/)) 



p{u,n) = n„ 



1 



u 



cot {^Kl{U)) 



^{U,n) = ^/2K-^nU 



2 -1772 ^[2 + U 
1 



2Kl{U) 
2 



2j2~U/2 T^[2 + U 



cot {^Kl{U)) 
2Kl{U) 



(105a) 
(105b) 

(105c) 

(105d) 



where the braid constants A„ and n„ are given in Table |2l These funct ions are plotted in Fig. [8] 
for the trefoil topology. We recall that a and /3 measure the first order perturbation to the 
internal force in the loop, see Eq. (51); p and (j) measure the infinitesimal rigid-body translation 



and rotation of the loop, respectively, see Eq. (53) 



At this point, we have expressed the internal parameters of all three regions, "^l and 
^B, as a function of the dimensionless loading U and knot type n. By plugging back these 
parameters into the solutions for the tail, loop and braid regions derived in Sections [5] [6] and |7] 
one defines a unique and smooth solution of the Kirchhoff equations representing a loose knot: 
we have eventually solved the problem formulated in Section [2j The solution, indexed by the 
dimensionless loading parameter U, is visualized in Fig. [9] 



8.3 Validation by direct numerical integration 

In order to check the analytical results, we have performed numerical simulations of knotted 
rods in the finite e case. Kirchhoff equations ([?]) were integrated numerically to find equilibrium 
configurations of a rod of finite thickness, knotted in an open trefoil, with a simplified contact 
topology (isolated contact points). Numerical continuation was then used to reduce the rod 
thickness and the leading orders for the position, tangent, internal moment and force were con- 
firmed, up to a small error due to the small penetration taking place in this approximate contact 
topology, see Appendix] A. 2 1 



8.4 Instability of the knot 

We have formulated the problem of finding the equilibria of the knotted rod as a set of linear 
equations expressing matching conditions between tail and braid, and between loop and braid. 
This linear system is regular except for some critical values of the loading U: 

|C7| = 2,V6,4, 
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Figure 8: Loop internal parameter s a, (3 , p and 
for a trefoil knot {n 
and 



as functions of reduced loading parameter U, 
a varies linearly with U . The other parameters, /3, p 
all diverge dX U — ±2, as denoted by the dashed lines (red). This divergence points to 



1). By Eq. (105a 



the helical instability undergone by the tails as the applied torque approaches the critical value 
TJ =±2. 
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Figure 9: 3D representation of the solution, for different values of the twist parameter: (a) 
U = -1.65, (b) U = -1.1, (c) U = -.55, (d) U = 0, (e) U = .55, (f) U = 1.1 and (g) U = 1.65. 
Rotation of the loop about the y axis is visible here, and takes place with the angle e0 which has 
been plotted as a function of U in Fig. [s] These 3D plots are based on the analytical solutions 
of the matched asymptotic expansion in each region, and are rendered here with e = .2. Note 
that the continuity of the solution across the different regions is only satisfied asymptotically for 
small e; in this rendering, e is non-zero and there is a slight mismatch at the junction between 
tails (red) and braid (blue), and between braid (blue) and loop (black). The same holds for the 
inextensibility condition, which is only approximately satisfied in the figure. 
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Figure 10: (a) Hat angle of a cinquefoil (5i) knot locked by friction on a Nitinol rod with 
radius h = .44 mm. No force is applied on the tails (T = 0) which are perfectly straight, (b) 
Datapoint obtained by repeating the experiments with various knot types (open symbols for 3i 
knot, filled symbols for 5i knots) and rod radii. In addition, the single datapoint shown by an 
empty circle is extracted from the work of Tong et al. ( 2003 1 . The two straight lines are the 
predictions of our theory, Eq. (106), with no adjustable parameter. 



The first and lowest value, \U\ — 2, comes from the matrix Mt(L^) expressing the response of 
the tail in Eq. (|99|. As explained in Section 5.3 the tails become unstable with respect to helical 
buckling when the applied twist reaches the critical value \U\ — 2. This explains the divergences 
observed in Fig. [sj and the large rotation of the loop in the first and last frames of Fig. [o] when U 
approaches ±2. We have confirmed this instability by direct numerical solutions of the Kirchhoff 



2008); it is analyzed in more details in a follow-up 



equations for dynamic rods (Bergou et al. 
paper. 

The other critical values given in the list above, namely \/6, 4, . . . come from Eq. (104). 
They correspond to the well-known Michell's instability of a twisted elastic ring, also known as 
Zajac instability, see Michell (1890). The lowest critical value that makes the loop unstable, 
\U\ = ^/6, is still larger than that for helical buckling \U\ = 2: in the case of infinite tails, the 
tails of the knot always buckles first. In the case of tails with a finite length, the threshold for 
helical buckling becomes larger than 2; for short enough tails, Michell's instability eventually 
sets in first. 



9 Experiments 



We present some validation experiments for the twistless case {U 
complement those reported previously by Audoly et al. ([20071) . 



= 0). These new experiments 
They were performed using 

naturally straight, superelastic wires made of Nitinol, an alloy of nickel and titanium, of radii 
in the range h = 0.17mm to 0.44mm, and of length 2 m. In Section [QA] we study the angle 
of the tails in a knot locked by friction, when no force is applied on the endpoints of the rod 
(T — 0). In Section 9.2 we confirm the existence of the two symmetric openings in the braid 



region predicted by the theory, and study them quantitatively. 



9.1 Hat angle 

For our first series of experiments, we use the geometry in Fig. |10[ A trefoil or cinquefoil knot 
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is tied on a Nitinol rod and its ends are gently released. If the knot has been formed with a 
small loop, its radius increases as the tails slide along each other in the braid region, until it 
reaches an equilibrium value. If the knot has been formed with a big enough loop, it stays in 
equilibrium when the rod is released. In either case, this leads to equilibrium configurations such 
as the one shown in Fig. [TO| i. No force is applied on the endpoints, T = 0, as friction in the braid 
region prevents the loop from further expanding. We are interested in the angle if, called the hat 
angle, made by the tails in the presence of frictional locking. This angle Lp has been measured 
in experiments with rods of various diameters, both for simple (trefoil) and double (cinquefoil) 
knots. These measurements are summarized by the symbols in Fig. [TOl 

The analytical method derived in this paper has been established in the frictionless case, 
when the knot is held by a tension force T 7^ 0. As we show now, it can easily be extended to 
configurations of the knot locked by friction. Let us first analyze in order of magnitude how the 
equilibrium radius of a locked knot depends on the coefficient of self- friction, which we call v. 
By our previous scalings, the internal force ns in the braid is of order 1/e, and the braid length 
is of order e. This implies that the contact force per unit length is of order By Coulomb's 

law, the tangential contact force per unit length is of order v/e^ , and the total friction force 
integrated along the braid is ^ v /e. The internal force is now zero in the tails: like the external 
tension T in the frictionless case, the integrated friction force must balance the internal stress 
in the loop to allow global equilibrium. Therefore v/e must be comparable to the loop stresses, 
which are of order 1. We conclude that v = 0{e). In other words if friction is weak, 1/ 1, 
the radii R compatible with equilibrium are those such that e = ^Jh/ R = 0{v); if friction is 
not weak, ly = 0(1), then e = 0(1) too, meaning that equilibrium configurations of the knot 
are tight and the present theory does not apply. In the experiments reported here, the friction 
coefficient was independently measured as « 0.1; this is consistent with the loop radii observed 
at equilibrium, which are such that .05 < e < .20. This reasoning shows that we must view the 
friction coefficient as a quantity of order e in our theory in order to consistently account for 
frictional locking. 

Knowing that v must be seen as a quantity of order e, it is now straightforward to adapt our 
matched asymptotic expansions. Indeed, the braid solution is not modified at dominant order by 
friction. The generic loop solution is obviously not modified cither. The only change concerns the 
tails whose loading geometry has changed; its ends are now free of any applied force or moment, 
and so both and niT are everywhere zero. As a result, U = and the perturbed tail solution 
given in Section [5] has to be replaced by perfectly straight tails. The functions xt{z) and y^iz) 



are affine functions of ^ and the expansion (43 1 is recovered, but with arbitrary coefficients Xt, 
Xip, Yt and Y^. When the pulling force T is zero, the axis z no longer plays a special role and the 
system becomes invariant by infinitesimal, rigid-body rotations about the y axis and translations 
along the y axis. We can use the rotation to make the tails perpendicular to the x axis; this 
amounts to set Xlp = by convention. Similarly, the translation can be used to set Yt = Q 
by a convenient choice of origin. With these conventions, Xt and Y^ can be chosen arbitrarily 
while X'j, and Yt are zero. This change can be accounted for by redefining the matrix Mt([/) 



in Eq. (|44| as follows: Mt = {{1,0}, {0,0}, {0,0}, {0, 1}}. This is the only change required to 
account for friction and self-locking. 

The matching procedure can then be repeated with the new tail matrix My. The hat angle 
is given by (yS = 2 |yy(0)| and we find: 

= n„ = V2n„e. (106) 



The value of n„ depends on the knot type and is given by Tab. |2l The prediction appears 



in Fig. 10 as the two straight lines for n = 1 and n = 2. There is a good agreement with 
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Figure 11: (a) Trefoil knot tied in a Nitinol rod of radius h — .44 mm and viewed from side. 
Loop radius is i? = 7.95 cm and e = ^/h/R = .075. (b) Close -up view of the same experiment 
revealing the two symmetric openings around the center of the braid. The bars of length Az ~ 
mm indicate the predicted apparent length of the openings, with no adjustable 
parameter, (c) Prediction for the apparent length of the openings is based on the fact that the 
endpoints Wa and a'^, are such that jy'' — y'^| = 2 /i. 



experiments for both knot types, especially in the range e < 0.1 — for larger values of e, the 
loose knot approximation appears to be less accurate, which is not surprising. 



9.2 Apparent length of openings in braid 



The second validation experiment concerns the two symmetric openings in the braid region, 
corresponding to (Te < | cr | < CTp in Eq. (88 1. The presence of these symmetric openings has been 
reported in our previous experiments, see Audoly et al. (20071. Here, we propose a quantitative 
validation: we consider the apparent length of these openings when the knot is viewed from 
the side, and compare the experimental measurements to the theoretical value. An experiment 
where these openings are visible is shown in Fig. |11| 

To predict the apparent length of the openings from our theory, we not e th at the endpoints 
of this region correspond to \y'' — 



2 h, as shown graphically in Fig. 
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In view of the 



rescahngs (33) and (76b I, this corresponds to v{a) = +1 (endpoints of the apparent opening on 
the positive side of the z-axis) or to v{a) = —1 (opening on the negative side). From Fig. [bJ), 
the function v{a) has a maximum slightly above 1 in the interval (fe < (t < (Tp. We call aa 
and the two roots of v{a) = 1 located on both sides of this maximum. Using the numerical 
solution of the universal difference problem given in Section 7.7.2| numerical root-finding yields 
the values of a a and a'^ for a trefoil knot, as well as their separation Act: 

Wa = 1.771, a'a = 2.230, Act = (ct^ - Wa) = .459 

In physical variables, this corresponds to an apparent [f]len gth of the openings Az = V2hR Act 
whose numerical value is Az = 3.86 mm in this particular experiment. This prediction is shown 

^^Notc that the actual length of the opening is much larger than the apparent length observed when looking 
along the x axis: in this particular experiment, the actual length of each opening is V2hR{ap~ae) = 19.5 mm. 
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by the two horizontal bars in Fig. |11| 3 and is in good agreement with the experiments, with no 
adjustable parameter. 



10 Conclusion 

We have considered the equilibrium of a knotted elastic rod under combined twist and tension. 
In general this problem should be expressed as a self-contact problem in 3D elasticity with finite 
strains and rotations. In this paper, we have considered the case where the theory of thin elastic 
rods is applicable, namely h <^ yj B /T where h is the small filament radius, T is the applied 
tension and B the bending stiffness. A crucial remark allowed us to derive analytical solutions of 
this problem: the assumption h <^ B /T warranting applicability of the thin rod model implies 
that the centerline is almost straight in the contact region. As a result, we could linearize the 
Kirchhoff equations in the region of contact, and formulate an equivalent contact problem with 
a fixed external obstacle. Our solution features a non-trivial topology of contact consisting of an 
interval flanked by two isolated points. 

We stress that, for all values of the parameters, the linearization of the equations in the 
region of contact is an approximation that is at least as good as the thin rod approximation 
itself. This remark could be applied to solve other geometries of rods in self-contact, such as the 
coiled configurations of elastic rings. This problem has been studied by numerical continuation 



by Coleman and Swigon ( 2000 1 . We expect that it can be solved by the same analytical method 
as the knot. One of the benefits of an analytical solution over a numerical one is that it captures 
the behavior of the equilibria for arbitrary values of the small thickness h and not just for specific 
values of h. 

Another interesting perspective opened up by the present work concerns the instability ob- 
tained for U = ±2, when helical bucking sets in in the tails. In a follow-up paper, we shall study 
this instability in details. Based on a refined version of the present theory, tailored to the case 
U « ±2, we shall show that the instability, which is driven by the tails, is strongly affected by 
the presence of loop; we also study what happens above the instability threshold. 



A Ruling out alternative contact-set topologies 



Contact problems are often solved by first inferring the topology of the contact set. Validation 
of this assumption requires checking that there is no penetration and that the contact pressure 



is everywhere positive. The approach we took in Section 7.7.2 is different as the topology of the 
contact set was found from constrained numerical minimization with no a ■priori assumption. 



We found an interval of contact fianked by two isolated points, see Eq. (88). Here, we investigate 



two alternative, simple contact topologies, namely a single interval of contact or three isolated 
points, and show that they lead to inconsistencies (negative pressure and/or self-penetration). 



A.l A single interval of contact 

Assume that the contact set is the interval a G [— CTi; CTi]. In this interval, the difference rod lies 
on the surface of the cylinder and can be parameterized as 



u{(j) = cos (/>(cr) , v{a) — sin <j){(7) . 



(107) 



Introduce the azimuthal vector = (— sin 0(tT), cos 0(a)). By deriving Eq. (1071 three times, 
we find 
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Now, the discontinuity of the the third derivatives {u"',v"') at the hft-ofF point is given by the 
point-like contact force, which is perpendicular to e^(ai) in the absence of friction. Therefore, 
{u"',v"') ■ 60 is continuous across cti, even though {u"',v"') is not. In addition, note that the 



asymptotic boundary conditions (82 1 for the braid imply that («"',«"') — (0,0) beyond the last 
contact point. We conclude that {u'" , v'") ■ is zero in the left neighborhood of the lift-off point, 
noted (7i^, which implies: 

0"'(ai-) = cj,'\wr). (108) 



By Eq. (78 1, the contact pressure can be found by deriving Eq. (107) four times. This yields 



p{a) = <f)'{a)^ - 3 0"(ct)2 -4(j)'{a) 0"'(ct). Combining with Eq. ( |108| , we compute the contact 
pressure at cti^: 

p(ffi-)--3 0''(ai-)-3 0"'(ar). 

This pressure is negative, which shows that the assumed topology of contact is inconsistent. 
Note that the pressure is negative in a region where the physical solution has openings, which is 
consistent. 

A. 2 Three isolated points 

Here, we assume that the contact set is composed of three isolated points: by symmetry, it is 
of the form V — {— cti} U {0} U for some cti > 0. As we shown now, this simple contact 

topology can be solved analytically and gives rise to residual penetration. We derive the solution 
on the positive part of the axis, ir > 0; the solution on the negative part can be found using the 
symmetry conditions ( [80| . 

Over the interval < a < +oo the contact pressure p is given by a Dirac function, noted 
Sd, centered at Wi: pia) — Pi5d(ct — cti). Noting — u{ai) Pi and Pi — v{ai) Pi the 



components of the contact force, we can write Eqs. (78 1 as u""{a) — a/2 P" <5z)(ct — ai) and 
v""{a) — \f2P^ (5d(<7 — (Ti). The general solution of these equations satisfying the asymptotic 
conditions (1821) reads 



= (Co + Ci ^) + ^/2 Pi e(ai - a) ^^^y^ (109a) 
v{a) = (^Co + ^ - y) + ^Pl e(ai - a) ^^l^L (109b) 
where Co, Cii Co ^"^d are constants of integration, and is the Heaviside function defined by 



0(a:) — for x < and 0(a;) — 1 for a; > 0. The expressions (109) are valid over the interval 
< fj < 4-00. Note that the right-hand sides are piecewise polynomial functions of a that are 

smooth; their third derivatives undergo a jump {\/2P^, \/2 P^) at a = ai. For a > ai the 
function Q is zero and u and v are given by the first terms in parentheses, while for < (f < cti, 
we have = 1 and u and v are given by third order polynomials. 

The seven unknowns of the problem (Co, Cij Co: Cii -P": -F'l'i'^i) can be found by solving the 
seven following equations: 

v(0) = 0, u'(0) = 0, v"{0) = 0, (110a) 
w(0) = 1, (110b) 
{P^,Pn-{u'{ai),v'{<Ji))^0, (110c) 
u'^(ai) + v'^{ai) = l, u(ai)u'{ai) + v(ai)v'(ai) ^0. (llOd) 



Eq. (110a I comes from the symmetry conditions near the center of the braid. Eq. (110b I comes 
from v{0) — and from the contact condition ^^(0) + ^^(0) = 1, which imply u(0) — ±1; 
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Figure 12: Radial distance w{a) = ■\/m^(o') + v'^{a) in the case of three isolated points of contact. 
The non-penetration condition w > 1 is violated near center. 



we consider u(0) = +1 only as the case u{0) = — 1 can be recovered by applying a symmetry 
X I— *■ (—x). Eq. (110c) warrants that the direction of the contact force is perpendicular to the 
tangent in the absence of friction. Eq. ( |110d ) expresses the fact that the rod has to be tangent 
with the cylinder at ai. 



In a first step, solve Eqns. (110a I and (110b I which are four linear equations for the variables 
Co, Ci, Co and P^;. This yields 



Co = l- 



3V2 



Ci = 



1 ^1 



PI = 



1 



(111a) 



Plugging these relations into Eq. ( 1 10c ) , we obtain a linear equation for Ci whose solution reads 

C;-ai-ai3(p-)2. (111b) 



Substituting into Eqns. ( llOd I, we find two polynomial equations for the two remaining unknowns 



and CTi, which have a unique real root 



o\ = 



(2 +77)3/4 

21/4 



2.661, 



P« = - 



21/6 ff^ 



5/3 ■ 



(Ulc) 



Eqns. (1091 and (111 I define in closed analytical form the unique braid solution having three 



isolated points of contact. 

Consider now the Taylor expansion of the distance function w = \Ju^ + near the center 
of the braid, w(ct) = 1 + ^w"(0)ct^ + . . . The coefficient w"(0) can be calculated as w"(0) — 

— \/ a/T — 2 / (4-\/6) ~ —.082 and is negative. This shows that there is some penetration]^ w <\ 
near a = 0, as confirmed in Fig. [T2j the solution with three points of contact is unphysical. 
Penetration takes place around the central point of contact; this points to the fact that the 
correct topology is obtained by replacing this point with an interval of contact. 

^■^Relative to the cylinder radius, penetration is by about 1 %. As a result, the unphysical solution with three 
points of contact happens to be a good approximation to the actual solution derived in Section |7.7| for a trefoil 
knot (to approximate a cinquefoil knot, one would need five points of contact). For instance, cti CTp = 2.681, 
Co = -.8729 Si A„=i = -.8776, and Ci = 2.0882 n„=i = 2.0887. 
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